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Abstract 

Reduced order modeling (ROM) seeks to make the modeling of aeroelastic behavior 
practical by reducing computation time for design codes. Deforming grids are often used in 
aeroelastic problems to account for the deformation of the structure. Proper Orthogonal 
Decomposition (POD/ROM) is a ROM technique that operates in an index-space for 
computations, not accounting for changes in grid dynamics, and must be modified to 
reflect grid deformation properly. As a POD/ROM is developed, fluid dynamics modes 
are created based on the index relationship between grid points. The modes are then 
used to recreate the full-order solution. When the relationship between the grid point 
locations and the index space changes, the created modes are no longer valid because 
the new grid dynamics are not captured accurately. To investigate and account for the 
effects of grid deformation on POD/ROM, a new algorithm is developed that incorporates 
modifications to the usual formulation. Evaluation of the new algorithm is accomplished 
through application to three fluid-structure models, each adding an increased level of grid 
dynamics. 

The first, an oscillating cylinder, incorporates an analytical potential flow solution 
as a full-order model. This model completely decouples the grid motion from the solution. 
Deforming grid POD/ROMs require more modes for an accurate solution than POD/ROM 
for rigidly moving grids. Degradation in POD/ROM accuracy in transforming the ana¬ 
lytical solution to reduced order space is investigated and found to be due solely to grid 
deformation. 

The second model problem, an oscillating panel in cross flow with discretized po¬ 
tential flow equations for the full-order solver, evaluates forced grid motion. This model 
couples the grid motion to the resulting fluid dynamics but does not couple the fluid dy¬ 
namics back into the structural model. A POD/ROM of a static grid with a transpiration 
boundary is compared to a POD/ROM of a deforming grid. Deforming grid POD/ROMs 
are found to require more modes than static grid POD/ROMs for similar accuracy levels 
for this model problem. In addition, for deforming grids, POD/ROMs are less accurate 


XVI 



when the grid deformation is significantly altered from the deformations seen in the devel¬ 
opment of the POD/ROM. As a result, a Multi-POD technique has been developed that 
evaluates the relative grid motion between how the POD/ROM was created and how it is 
executed. The Multi-POD technique determines the current relative grid deformation and 
selects the best POD/ROM from those available. 

Finally, Multi-POD is applied to a pitching and plunging airfoil with both forced 
and free dynamics. In the free dynamic case, the model fully links both structural (grid) 
dynamics and fluid dynamics. POD/ROMs are trained with forced grid deformation. In 
cases with free grid deformation, the Multi-POD technique is able to reproduce the free 
deformation solution very accurately, switching between POD/ROM as necessary. 


xvii 



TECHNIQUES FOR REDUCED ORDER 
MODELING OF AEROELASTIC STRUCTURES 
WITH DEFORMING GRIDS 

I. Introduction 

In the design of aerospace systems, designers and operators are looking for ways to 
computationally model prospective vehicle configurations quickly. Designers typically use 
low-fidelity models to generate vehicle designs rapidly, allowing the designer to evaluate 
a large number of design variations. As the design is finalized, the designer moves to 
more complete models. However, as increasing performance objectives drive designs to 
become more complex, the designer is faced with a dilemma: low-fidelity models cannot 
accurately model even the initial design variations. Higher fidelity models are desired, 
but are impractical for design due to their large computational expense in both time 
and hardware. To model a configuration fully, time integration with Computational Fluid 
Dynamics (CFD) is prohibitively expensive owing to the large number of degrees of freedom 
(dof) in the problem. Reduced-order modeling (ROM) is one means of developing higher 
fidelity models that are more efficient to solve, and therefore more attractive to designers. 
By reducing the dof (independent variables) of the problem, computational times can be 
decreased by orders of magnitude while still maintaining solution accuracy. 

Modeling has two basic attributes: fidelity and complexity. Fidelity is a model’s 
ability to predict the behavior of the system correctly. Sometimes called accuracy, it is 
a reflection of the degree to which the physics of the modeled system are described by 
the mathematics contained in the model. Complexity is how difficult the model is to 
generate and solve, often represented by the number of dof in the discretized, constituent 
equations. With the growing complexity of computational methods, model reduction has 
been studied for the analysis of very large problems. For problems with complicated 
physics, models with higher fidelity usually have a higher complexity. Model reduction 
seeks to produce an optimal trade-off between fidelity and complexity, producing a suitably 
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accurate model, while minimizing computational time. One such technique is Proper 
Orthogonal Decomposition (POD). 

POD was introduced by multiple sources in the mid 1940’s (1, 2). The concept has 
been applied in a wide variety of fields: fluid mechanics, uncertainty analysis, image pro¬ 
cessing, signal analysis, data compression, process identification and control in chemical 
engineering, and oceanography (3, 4, 5, 6, 7, 8, 9, 10, 11, 12). In its early development, POD 
was presented as a mathematical method of correlating statistical data (13, 14). Lumley 
(15) introduced the idea of applying POD to turbulent fluid flow in 1967 as a modeling 
technique. Reduced-Order Modeling via Proper Orthogonal Decomposition (POD/ROM) 
for turbulent modeling is still in use today (16). In its application for reduced order models, 
POD/ROM was based on earlier forms of eigenmode analysis. Using eigenmode analysis as 
a means of generating a reduced order model for fluid dynamics first appeared in the late 
I980’s and early I990’s. Dowell presented several papers on the development of eigenmode 
analysis for ROMs, and his work (17, 18, 19) represented a literature review of the subject 
area. Eigenmode analysis for ROM of fluid problems has been applied to several areas of 
interest: classic incompressible and potential flows, compressible potential flow, compress¬ 
ible Euler flow, and unsteady viscous flow. One of the initial papers of eigenmode analysis 
of classic incompressible and potential flows was presented by Hall (11) in 1994. Hall de¬ 
veloped eigenmode based reduced order models for aerodynamic representation of inviscid, 
incompressible vortex-lattice models, to represent two and three dimensional flows about 
airfoils. Mahajan, Blahle, and Dowell (20) calculated the eigenvalues of an aeroelastic sys¬ 
tem modeled with the full-potential equation to investigate aeroelastic stability. Building 
on the above works, Romanowski and Dowell (21) generated the first successful generation 
of ROM’s for linearized Euler flows and showed a definite improvement in calculation time 
for complex aeroelastic problems. Florea and Hall (22) then moved on to create reduced 
order models in the time domain for linearized potential flow about isolated airfoils. In 
1999, Hall, Florea and Dowell employed the Lanczos procedure (23) to create a reduced 
order aerodynamic model in the frequency domain for linearized unsteady potential flows 
in turbomachinery. POD/ROM was used in 2001 by Thomas, Dowell and Hall for model 
problems in three dimensions (24). 
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POD/ROM is a spectral method using modes created by taking snapshots (time 
dependant data samples from CFD program or experiment, captured at increments in 
time and representing the dynamics of the fluid system) of the unsteady flow-field data. 
These snapshots are used to generate fluid modes using an eigen-problem method. If the 
snapshots sufficiently represent the fluid dynamics of the problem, the POD/ROM will 
be accurate. The results of POD/ROM could then be used to analyze similar, but not 
identical, conditions to the initial model. The range of accuracy around the original model 
represents the fidelity of the POD/ROM. 

POD/ROM was first applied to aeroelastic problems in 1996 by Romanowski (25) 
as a means of generating the basis set for the linearized Euler equations. Dowell, Hall 
and Romanowski (19) used this implementation in the evaluation of a plunging and ro¬ 
tating airfoil. Later, Ly and Tran (26) used POD/ROM as a method of generating a 
Galerkin reduced order model for a chemical vapor deposition reactor, showing an effi¬ 
cient approximation of solution for compressible viscous flows coupled with energy and 
species equations. lollo, Lanteri and Desideri also developed Galerkin approaches for solv¬ 
ing viscous flows (27, 28). In 1999, Beran, Huttsell, Buxton, Noll, and Osswald (4) used 
POD/ROM to generate a reduced order model for the calculation of limit cycle oscillation 
of a convection-diffusion-reaction equation. In 2000 and 2001, Beran, Pettit and Mortara 
used POD/ROM for nonlinear panel response (29, 30). POD/ROM was also expanded to 
coupled fluids and flow control problems by Arian, Fahl, and Sachs (31). 

In 2000, LeGresley and Alonso (32) used POD/ROM for the optimization of an airfoil 
design. This effort used mesh deformation with POD/ROM, though it did not specifically 
address it. An Euler solver was used to generate steady state snapshots for various airfoil 
shapes. A POD/ROM was developed from these snapshots and a least squares optimization 
used to change the airfoil shape based on input parameters to the POD/ROM. LeGres¬ 
ley and Alonso were able to change their airfoil shapes to meet new input parameters 
when the airfoil shape changes were small. However, they did not examine any aspect of 
the deforming grid, nor how it affected the accuracy and convergence properties of their 
POD/ROM. 
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1.1 Research Objective 

The POD/ROM research discussed used either small deformation theory or rigid 
grids. It is the objective of this research to demonstrate how POD/ROM of fluid dynamics 
problems is affected in terms of efficiency and effectiveness by the use of deforming grids. 
Grid deformation is used in CFD to model a variety of moving boundary problems, such 
as dynamic interaction between airfoils and flaps, and deforming aeroelastic structures. 
POD/ROM computations are conducted in computational (or index) space and do not 
directly account for grid deformation. Snapshots taken of the flow-field data are arranged 
in column (index) form. The order of the data in the column is irrelevant to the process, 
but must be consistent. In any problem with a deforming grid, the physical location of 
any data point changes from one snapshot to the next, even though its location in the data 
column does not. This makes the deforming grid POD/ROM more complex than that of 
a static grid POD/ROM where the relationship between index space and the physical grid 
locations is fixed. 

To examine this effect a model problem with a steady analytical solution is evaluated: 
the analytical solution of potential flow about a uniformly translating 2-D cylinder in still 
fluid (see Appendix A. for full model development). An analytical solution decouples the 
grid motion from the fluid dynamics. Two grid cases are selected. The first case is a rigid 
grid fixed to the cylinder and translating with it (no grid deformation). The second case has 
a fixed outer boundary and the grid is attached to the cylinder as it is translated and rotated 
within the boundary, causing grid deformation. For each grid case, POD/ROM is applied 
to snapshots of data obtained by imposing the analytical solution on the grids at selected 
increments in time. With the analytical solution decoupling the grid dynamics from the 
fluid dynamics, any difference in the POD/ROM of the two cases can be solely attributed 
to the grid deformation. The resulting modes created from the snapshots represent the 
differences in grid dynamics for the two cases. The relative contribution of each mode can 
be estimated by its percentage contribution to the total eigenvalue energy. This gives an 
idea of how important individual modes are to the overall POD/ROM solution. 

The resulting POD/ROM of the two cases are significantly different. By plotting 
the modal contribution verses the number of modes, the relative importance of each mode 
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# Modes 

Figure 1.1 POD Eigenvalue Magnitude 

can be determined. Differences in the slope of modal contribution verses number of modes 
highlight differences in POD/ROMs. If modal contribution drops quickly, only the initial 
modes are of significance. If the slope is more shallow, more modes are necessary for an 
accurate solution. The POD/ROM for the rigid grid, POD/ROM/RG, has a first mode 
with 100% of modal contribution (100% of the total energy of the eigenvalues) and the 
remaining modes are constrained by machine precision. This indicates that the ROM of 
the analytical solution on a rigid grid is exact. The POD/ROM for the deforming grid, 
POD/ROM/DG, has a first mode with 99.7% modal contribution. The remaining modes 
decrease in energy quickly, but are still significant compared to the non-dominant modes of 
the POD/ROM/RG (Fig. 1.1). Therefore, the POD/ROM/DG is not an exact ROM of the 
analytical solution. Additional modes are necessary for an accurate solution. By plotting 
contours of the stream function variable, T, the analytical solution can be compared to the 
POD/ROM/RG solution. The POD/ROM/RG with one mode is exact (Fig. 1.2) when 
compared to the analytical solution. The single mode necessary for the POD/ROM/DG is 
a scalar percentage of the exact solution. The mode can be captured at any time increment 
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(a) Analytical Solution (b) POD/ROM with 1 Mode 

Figure 1.2 Contours of Potential Variable (Rigidly Attached Grid) 

with a single snapshot. In contrast to the POD/ROM/RG, the POD/ROM/DG requires 
at least 7 modes for a 99% accurate solution. If fewer modes are used the POD/ROM/DG 
solution degrades. For contours of ll/, a 4 mode POD/ROM/DG compares poorly to the 
analytical solution (Fig. 1.3). The principle mode is incorrectly mapped to the deformed 
grid and results in a skewed solution. The additional modes necessary for an accurate 
solution map the principle mode to the deformed grid. The additional modes necessary 
for a POD/ROM/DG over a POD/ROM/RG are solely caused by the addition of grid 
deformation. 

To use a POD/ROM successfully in cases where the deforming grid interacts with 
the fluid dynamics, the magnitude of the increase in the error needs to be understood 
and quantified. The model problem described has the fluid dynamics decoupled from the 
grid motion, however, it clearly demonstrates that grid deformation affects the accuracy of 
POD/ROM in an analytical problem. In more realistic cases, the fluid dynamics and grid 
dynamics are coupled, such as in a pitching and plunging airfoil. The coupled nature of 
the dynamics precludes a simple method of evaluating the error due to the deforming grid. 
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(a) Analytical Solution (b) POD/ROM With 4 Modes 

Figure 1.3 Contours of Potential Variable (Deforming Grid) 

However, the error resulting from a POD/ROM/DG may be estimated by comparing the 
results to a case for which the grid is static. 

POD/ROM are typically accurate about some standard state. Parameters are varied 
about a standard state, such as angle of attack for an airfoil, and provide the POD/ROM/DG’s 
range of accuracy. In a POD/ROM/DG, the grid deformation itself becomes a parameter. 

As the POD/ROM/DG is developed by collecting snapshots (called the training period), 
the modes created are based on the different grid states seen in those snapshots. To be 
robust, POD/ROM/DG must have some range of accuracy over variations in grid states 
seen in the snapshots during training. This adds an additional complexity to the creation 
of a POD/ROM/DG. If the POD/ROM/DG is applied to a grid state not included in the 
training set, regardless of whether the rest of the parameters are identical, the accuracy is 
affected. The magnitude of the allowable change in deformation defines the robustness of 
the POD/ROM/DG. 

1.1.1 Thesis. The thesis of this dissertation is that POD/ROM can he used to 
simulate fluid flows accurately on deforming grids. 
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1.1.2 Approach. The research was decomposed into four steps. The first step 
was to evaluate POD/ROM with deforming grids, where the fluid solution was independent 
of the grid. Potential flow around a translating cylinder was selected because it had an 
analytical solution. This allowed a complete separation between fluid and grid dynamics. 

The second step was to develop a ROM for the unsteady full-potential solver. The 
solver used was the Shankar (33) implicit scheme. The scheme was relatively simple code, 
allowing more attention to be focused on the deforming grid and its application to a 
POD/ROM. The POD/ROM implementation was based on the RAPOD computer pro¬ 
gram developed by Reran (4). RAPOD was modified for a fully implicit scheme based on 
the work of Dowell, Hall, etc. (34). The Shankar solver was validated using both steady 
and unsteady experimental data. 

The third step was to analyze the effects of the deforming grid on the POD/ROM 
solution, where the grid dynamics would affect the fluid dynamics. The oscillating panel in 
cross flow was selected as the model problem. Analysis included comparison to a static grid 
POD/ROM as well as full-order solutions based on the full-potential solver. Techniques 
were explored to evaluate the robustness of POD/ROM based on changes in the grid 
deformation. 

The fourth step was to examine a fully coupled structural and fluid dynamics prob¬ 
lem. The pitching and plunging airfoil was chosen as a model problem and the techniques 
developed from the oscillating panel were applied to determine their effectiveness in this 
more general case. Finally, POD/ROM was tested for a practical problem of free pitch 
and plunge to determine its utility with deforming grid aeroelastic problems. 

1.1.3 Research Questions. Several research questions are addressed in this work: 

1. What decrease in accuracy can be expected solely due to the grid deformation for a 

POD/ROM? 

2. What metric can be used to evaluate the accuracy of a POD/ROM on a deforming 

grid? 
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3. Can a POD/ROM constructed for a deforming grid be used for a grid with different 
deformation properties? 

1.2 Contribution 

The contribution of this research to the body of knowledge is in expanding the 
application of POD/ROM to a wider variety of problems via relevant model problems. 
This has been accomplished using a deforming grid that allows for the analysis of moving 
airfoils, flaps, and structures. 

1.3 Document Organization 

The remainder of this document is organized according to the approach described in 
Section 1.1.2. Chapter Two examines the application of Proper Orthogonal Decomposition 
to the CFD model chosen. Chapter Three looks at deforming grids and error estimation 
for POD/ROM. Chapters Four and Five address two application problems: the oscillating 
panel and the pitching and plunging airfoil. Chapter Six summarizes the research findings 
and provides ideas for future work. The appendices detail the translating cylinder model 
problem, the theoretical background on POD/ROM, and the derivation and validation of 
the flow solver. 
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II. Proper Orthogonal Decomposition 

Proper Orthogonal Decomposition (POD), also known as Karhunen-Loeve decomposition 
(1), is an empirical spectral method. Spectral methods use basis to approximate the full- 
order solution. The basis can be developed through a variety of techniques: trigonometric 
functions, polynomials, or wavelets for example. The basis is used to generate constant 
modes that are constant in space and are multiplied by coefficients that are time-dependent 
to create an expansion to the time-dependent, full-order solution. That is, if cn is a column 
array representing a set of N variables, it can be represented by a set of M modes, 0^, 
and a set of reduced-order variables, Cjm {t), such that 

M 

^ (^) ~ 'y ^ (^) 4’m- (2-1) 

m=l 

where N » M. It is the development of the modes, 0^, that differentiates one spectral 
method from another. 

In POD, the basis functions are created from a set of observed data, herein called 
snapshots. As a result, the POD basis can be made optimal in terms of minimizing the 
error in recreating the observed data (see Appendix B). This does not mean it will be the 
optimal basis for all, or even most, full-order cases, only for the snapshots from which the 
modes are developed. This is one of POD’s weaknesses: to be accurate the snapshots must 
fully represent the dynamics of interest in the fluid model. For instance, data from a CFD 
model would need to reflect the complexity of the model problem and be complete enough 
to capture all relevant fluid structures, including non-linear effects. The creation of the 
snapshot is called the training period. The training period comes at a computational cost 
and the more complete the set of snapshots, the more expensive the training. However, if 
the snapshots are not complete the ROM will not be robust to changes in parameters that 
govern the system behavior. As Beran and Silva (35) noted in their development of POD, 
the trade-off can be favorable when, after the initial computational investment, a compact 
ROM can be constructed which can be used many times, and is valid over a useful range 
of system states. For deforming grid cases, the issue lies in creating a POD/ROM that is 
useful over a range of deformations. 
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In POD, the snapshots are linearly combined to form a smaller number of basis 
vectors. That is, for K snapshots of uo, 

K 

4’m ~ 'y ( 2 - 2 ) 

k=l 


where is the contribution of the kth snapshot to the mth basis vector. In matrix form 
Eqn. 2.2 becomes, 

$ = SV, (2.3) 


where, 


and 


T T T 

01 02 0m ’ 
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(2.4) 


(2.5) 
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( 2 . 6 ) 


is known as the modal matrix, an {N x M) matrix containing the functions 0^. The 
snapshot matrix, S, is an {N x M) matrix containing M snapshots of u (t), with the data 
stored as column vectors in the matrix. The transformation matrix, V, is the (M x M) 
matrix that maximizes the projection of the snapshot matrix onto the POD basis. V must 
be properly scaled to be orthonormal. This leads to the eigenproblem. 


S'^SV = VA, 


(2.7) 


for eigenvectors, V, and eigenvalues, A = diag (Am) (see Appendix B. for the full devel¬ 
opment of V). In practice, fewer modes are used than the total number of snapshots. 
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M < K. These are selected based on the relative size of the modal eigenvalues, Am, in a 
technique called truncating. Each eigenmode is normalized and compared to the sum total 
of all eigenmodes for its modal contribution. 


2.1 Full-Potential Equation 

The scope of the current work uses the full-potential equation for the governing 
equation (see Appendix C. for full development). The full-potential equation is an approx¬ 
imation to the Euler equations for which one assumes no rotation in the flow, and thus no 
entropy production (36). The potential equation has one fluid variable, T, defined as. 


dx ’ 
dy ■ 


( 2 . 8 ) 

(2.9) 


where u and v are fluid velocities. In a 2-D, body-fitted coordinate system represented by 
T = t, ^ ^ {x,y,t), and y = r]{x, y, t), the unsteady full-potential equation is written as 


P 

J 


T 


+ 




( 2 . 10 ) 


where the (i?)^ refers to the derivative of R with respect to U and V are the contravarient 
velocities, J is the Jacobian of the metrics and the density, p, is represented as 


P = 


= 1 
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Ml [2J/, + {UF Tj + (E + J/, - 1] 


1 
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( 2 . 11 ) 


The above equations can be solved by a variety of schemes. The scheme selected was 
developed by Shankar, et, al. (33) and has the benefit of being very stable and easy to 
implement with POD. It uses an implicit Newton’s method to solve Eqn. 2.10 for a robust 
and efficient numerical solver with effective treatment of boundary conditions. In summary, 
the scheme is as follows; Eqn. 2.10 is represented as 


A(vk) = 0, 


( 2 . 12 ) 
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where is the unknown to be solved for at every grid point in the domain for time step 
(n + 1). The Newton iteration for the solution to Eqn. 2.12 is, 

(vk-T,) = 0, (2.13) 

where is the current guess for T at the {n + 1) time step. At convergence, AT = T — 
will approach zero. The Jacobian matrix, (§§), is large (N x N) but sparse and can be 
simplified using an approximate factorization creating two linear operators. 




L 


■n 


_ A rr ^ 1 P ^ 

l + Arl/|- + i^a22|- 

orj [3 J orj 


(2.14) 

(2.15) 


where jS = () > oii = 022 = ’lx + Vy- ^-i^d are (NxN) tridiagonal 

matrices. Eqn. 2.13 is then solved successive ^ and rj sweeps 


L^LyAT = R, 

(2.16) 

L^AT = R, 

(2.17) 

LyAT = AT, 

(2.18) 


where AT is an intermediate solution and R = F (T^.) from Eqn. 2.13. The above scheme 
has been validated for subsonic and transonic flows. It has been compared to steady and 
unsteady experimental data as well as other CFD codes (see Appendix D.). 


2.2 POD Implementation with Potential Equation 

POD can be implemented in a fully implicit manner in the above scheme. Following 
the work of Hall, et al. (34), the reduced order variable in Eqn. 2.1 becomes e, for the ^ 
sweep, 

AT = Te. (2.19) 

Substitution into Eqn. 2.17 yields, 

L^Te = R. (2.20) 
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Applying the transpose of $ to both sides the above equation, 


( 2 . 21 ) 

produces a square (M x M) matrix, where M is the number of modes in the 

POD/ROM. Because M « N, the matrix can be easily inverted to solve the ^ 

sweep, 


e = (2.22) 

^ (2.23) 

Tridiagonal techniques require a reordering from row major to column major between the ^ 
and r] sweeps. Likewise, Aik must be reordered before the next POD step. After reordering, 
the same POD steps are applied for the rj sweep; 


L^AT 

= AT, 

(2.24) 


= AT, 

(2.25) 


II 

> 

(2.26) 

e 

= ^ T^AT, 

(2.27) 


This can then be converted back to the full-order solution used in the Newton’s method, 

AT = T^AT. (2.28) 

Computational savings come from the inversion of and being less com¬ 

putationally expensive than their corresponding tridiagonal solutions of Eqns. 2.17 and 
2.18. The inversion of and could be done with a weighted pseudo¬ 

inverse for more accuracy. Weighting assigns more importance to modes that dominate 
the fluid dynamics. A weighted POD/ROM technique was presented by Willcox and 
Perairet in 2001 (37), that demonstrated balanced modes for a small deformation pitch¬ 
ing and plunging airfoil, showing good results. However, weighting techniques assume a 
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level of knowledge about the nature of the model problem and the way the modes affect 
the fluid dynamics. A user must know which modes to weight heavily. In addition, the 
weighting would require an additional post processing matrix multiplication, decreasing 
the efficiency. In the truncation method described, the weighting is binary, and based on 
the contribution of the eigenmodes. 

POD is an excellent technique for developing an optimal reduced-order model, from 
observed data. In practice, most of the computations for POD are done during the training 
process and do not impact the execution of the POD/ROM itself. This makes POD/ROM 
a computationally attractive technique as well. POD/ROM can be conveniently incorpo¬ 
rated to the full-potential equation. The application of POD/ROM with the full-potential 
equation is explored in later chapters through model problems. 
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III. Deforming Grid Error Estimation 

In most aeroelastic problems, the motion of the body is tightly coupled to the resulting fluid 
dynamics. As a result, it is difficult to separate the effects of grid motion from the coupled 
POD/ROM. A rigid grid, attached to the body in motion and moving with it, eliminates 
deformation. The resulting POD/ROM/RG can then be compared to a POD/ROM/DG 
to determine the exact effects of deformation. However, rigid grid motion is not possible 
in all problems, such as in the case of a moving flap attached to an airfoil. If an analytical 
solution is applied to a deforming grid, then the fluid dynamics can be decoupled from the 
grid motion and a POD/ROM created from only the grid deformation. However, analytical 
solutions are not available for many problems. To evaluate the effects of grid deformation 
on POD/ROM, some method of evaluating the quantitative effect of deforming grids on 
POD/ROM accuracy needs to be developed. Two issues dominate: the degradation of 
accuracy due to the deforming grid, and the robustness of POD/ROM/DGs. 

3.1 Qualitative Error Estimate 

To estimate degradation of accuracy, one technique of estimating the error for a de¬ 
forming grid POD/ROM is to limit the amount of deformation for a simplified problem. 
For instance, in an oscillating panel if the deflection is small enough, the number of de¬ 
forming grid points can be varied from a limited number near the panel, to more points 
further out (Fig. 3.1). This allows POD/ROM to be created with nearly identical fluid 
dynamics but differing deforming regions. The resulting solutions can be compared to 
their respective full-order solutions. The percentage of deformation and the resulting error 
for a fixed number of modes can give an estimate of the error resulting solely from grid 
deformation. 

3.2 Deforming Grid Metrie 

To examine the robustness of the POD/ROM on differing grids, a separate technique 
is needed. In examining the POD process, it is noted that the snapshots taken of the fluid 
dynamics during the training process represent the range of the modes created. That is. 


3-1 



only fluid structures captured in the training process, or their linear combinations, can be 
recreated by POD/ROM. Also, the snapshots of the fluid dynamics are arranged in index- 
space and therefore relationships between the respective snapshot of the grids are lost. 
Each grid is assumed to be identical to the previous. This is not the case with a deforming 
grid. Therefore, it is necessary to track the relative deformation of the grids during the 
training process. The training grids represent the range of grid motion. The grids captured 
during the training process, or their linear combinations, should have POD/ROM/DG 
with smaller grid deformation error. If the POD/ROM/DG is run on grids which differ 
significantly from the training grids, then the grid deformation error increases. 

To evaluate how closely the run grids match the training grids, a grid modal matrix 
can be created from the grid snapshots. Snapshot matrices of the x and y locations 
for the deforming grid are stored during the training process. The snapshots are based 
on the perturbation from a default grid to highlight the relative deformation. The grid 
snapshot matrices are used to create a grid modal matrix, in a fashion similar to the 
creation of the modal matrix of the POD/ROM/DG. The grid modal matrix represents 
the range of grid deflections during training. The grid modal matrix is correlated to the 
POD/ROM/DG, in that the grid snapshots are taken at the same points in time as the 
fluid snapshots. 

An error term can be generated from the grid modal matrix through a pseudo¬ 
inverse technique. Assuming that the grid can be represented by a linear combination of 
the training grids, 

X ss (3.1) 

where x is the exact value of the run grid. To compute x directly, a pseudo-inverse of 
is needed. Using a Generalized Inverse (38), and multiplying both sides of Eqn. 3.1 by the 
transformed matrix, produces a square matrix, on the RHS, 

^Ix = (3.2) 

that can be inverted. 

X = ^Ix. (3.3) 
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Finally, the error of the difference between the training grids and the run grid is determined 
by substituting x into the original approximation for x and subtracting it from the exact 
run grid x value, 

erroTx = x — (^g^g) ^ ^g^^ (3.4) 

The term $g (^g^g) ^ ^g represents how closely the run grid is represented by the 
training grid. The same technique can be used for the other spatial dimensions. 

Both the grid modal matrix and the pseudo-inverse need only be generated once; 
therefore the error evaluation does not add significantly to the overall computation time. 
Several metrics are considered for correlation to the POD/ROM/DG error: maximum 
percentage grid error, L 2 norm of the error, and the grid error scaled to the non-dimensional 
length unit of the problem. 

The grid error can be correlated to the error of the POD/ROM by means of running 
a test case. A POD/ROM is created at one set of grid deformations and then run at a 
second set of grid deformations. The error of the POD/ROM to the full solution for the 
second set of grid deformations is due to the misalignment of the modes to the grid point 
locations. By comparing the relative grid error and the POD/ROM error, the correlation 
between the errors can be determined. An acceptable level of POD/ROM error is then 
selected and the corresponding threshold of grid error established. Once the grid error is 
correlated to the POD/ROM error, the grid error is used as a threshold indicator for the 
accuracy of the POD/ROM. When the grid error exceeds a certain threshold, the accuracy 
of the POD/ROM has also exceeded the preset threshold indicating to the user that the 
POD/ROM should be replaced. 

These techniques are the basis for comparing deforming and static grid POD/ROMs. 
The deforming grid metric is used in the subsequent applications to evaluate robustness 
of deforming grid POD/ROMs and to determine the best available POD/ROM for any 
particular grid deformation. 
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(a) 100 percent Doman 


(b) 33 percent Domain 


Figure 3.1 Deforming Grid Domains 
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IV. Oscillating Panel 

To evaluate the issues of POD/ROM with a deforming grid, the nonlinear, discrete, full- 
potential equation in the simulation of unsteady flow over an oscillating panel was analyzed. 
This application coupled the fluid dynamics and the grid deformation. It required a re¬ 
duced order model to determine the effects of grid deformation, unlike the translating 
cylinder model problem (see Appendix A.) that evaluated only the accuracy of the ana¬ 
lytical solution projected into a subspace. The panel oscillation was forced and the fluid 
dynamics were not allowed to alter the grid motion. As a result, the predefined motion 
of the grid drove the fluid dynamics. The coupled nature of the fluid and grid dynamics 
made it impossible to evaluate modes solely caused by the grid motion. Therefore, the 
oscillating panel was modeled both with a deforming grid and a static grid transpiration 
boundary condition (TBC). This allowed a direct comparison of a POD/ROM on a static 
grids and deforming grids for a nearly identical problem, showing the effects of grid defor¬ 
mation on POD/ROMs. To evaluate accuracy of the POD/ROM/DG, the number of grid 
points that were allowed to deform was varied. Some types of problems preclude the use 
of the transpiration boundary condition or a rigidly moving grid. Therefore, varying the 
number of deforming grid points was used to provide a qualitative method of evaluating 
the accuracy of POD/ROM/DG. To evaluate the robustness of POD/ROM on deforming 
grids, POD/ROM/DGs were used on grids other than their training grids. The training 
grids were compared to the grids at which POD/ROM/DGs were run. The relative grid 
error allowed a measure of merit of the training process to be developed and correlated to 
errors of POD/ROM/DGs. 

4-1 Application 

Following the work of Pettit and Reran (39), the application of inviscid 2-D flow over 
an oscillating panel was used to examine POD with deforming grids. The flow-field was 
assumed to occur in an open channel above an infinite, segmented panel that nominally is 
in the y = 0 coordinate plane, except for a segment between x = 2.5 and x = 3.5 for which 
the panel shape was a smoothed parabolic defined by y (x, t). The chord length (measured 
from start to the end of the panel), c, of the panel was normalized to a value of 1.0. The 
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Figure 4.1 Model Problem 


channel was 2.25 chord lengths high and the inlet and outlet were 3.0 chord lengths from 
the center point of the panel. The domain dimensions were sufficient to provide a good 
solution for the freestream Mach number of 0.5. The Mach number was selected to keep 
the flow over the deflected panel subsonic to prevent shockwaves. A 80 x 30 grid was used 
for the model problem, with 20 grid points evenly spaced on the panel surface (Fig. 4.1). 
The surface was based on a smoothly varying function used for 2-D experimental wind 
tunnel tests (40): 

X (t) = cos (6) — [cos (6) — cos (30)], (4.1) 

Y (t) = sin (0) — [3 sin (0) — sin (30)] , (4.2) 


where 0 ranges from 0 < 0 < vr and r (t) was the time dependant amplitude of the panel 
at the midpoint, scaled to the chord length. The amplitude was varied based on a cosine 
function to give a zero panel velocity at the peak and zero deflections. 


r{t)=A 






(4.3) 


A was the maximum amplitude of the midpoint on the panel, / was the non-dimensional 
frequency of oscillation, t was the time, and tmax the maximum time for the test. 

The density and velocity were non-dimensionalized using the far-field density, p^, 
and velocity, Uoo- Pressure was non-dimensionalized using the far-field dynamic pressure, 
Poo (subsequent references to variables assume non-dimensional forms). 
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(a) Transpiration Boundary Condition Static (b) Deforming Grid 

Grid 

Figure 4.2 Computational Grids 

In the transpiration boundary condition case, the grid was static and undeflected. 
In the deforming grid case, the grid was smoothly deformed above the panel (Fig. 4.2). 
The panel was not allowed to dip into negative deflections. At the zero deflection point, 
the transpiration and the deforming grid cases had identical grids. 


Jf^A.l Boundary Conditions. Two boundary conditions were used to model the 
deflecting panel; the transpiration boundary condition and a deforming grid. The transpi¬ 
ration boundary condition is a small disturbance theory model that uses injected fluid 
(imparted velocity) through an undeflected transpiring panel to simulate a deflection. 
This allows the transpiration boundary condition grid to remain static. The transpiration 
boundary condition enforces the exact condition of impermeability of a deflected panel 
surface by the injection of velocity at the undeflected surface. 


—u- 


djh 

dx 


+ v = 


dys 

dt 


(4.4) 


^ is the slope of the deflected panel and ^ the time rate of change of the panel (41). 
This transfer of boundary conditions is identical to that employed in small-disturbance 
theory, where one assumes the regularity of the computed solution and small deformation. 
For the full-potential equation the boundary condition becomes. 


V ={u- Xr) (Ay)^ J, 


(4.5) 
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where Ay is difference in the y ordinates of the assumed deflected panel and the actual 
undeflected surface. The transpiration boundary condition assumes a quasi-static change 
in the grid, and therefore is limited to low oscillation frequencies for the model problem. 

The deforming grid boundary condition physically alters the grid between each iter¬ 
ation to match the new deflection of the panel. This requires that time dependant grid 
metrics be properly calculated and included in the governing equations. The deforming 
grid boundary condition enforces the exact condition of impermeability at the deflected 
panel surface, V = 0 (see Appendix C for full development). 

4-1-2 Full-Order Results- Full-order results are created during the training pro¬ 
cess to compare with the POD/ROM results. The model problem is non-linear in both 
maximum amplitude and oscillation frequency. Fig. 4.3 shows the periodic minimum 
normal force coefficient, Cn, after the flow is fully developed (approximately ten oscilla¬ 
tions) for a variety of oscillation frequencies and amplitudes. Three maximum amplitudes 
{A = 0.05,0.1.0.15) and nine oscillation frequencies are used (0.1 to 1.0), for 36 different 
combinations. The normal force coefficient is defined as, 

y- (Poo-P) 

Cn = ■ (4.6) 

Q.OO 

As the amplitude is increased, the problem moves from subsonic to transonic flow (at about 
A = 0.17). This study is interested in subsonic flow fields and so the amplitude is limited 
to 0.15. The increasing slope of Cat verses frequency shows increasing non-linear behavior 
with maximum amplitude. Breakpoints in the slope at frequencies of 0.4 and 0.65 show 
changes in linear behavior with frequency. To be accurate, a POD/ROM must be able to 
capture these non-linear behaviors. 

In addition matching derived variables, POD/ROMs must be able to simulate the 
full-order solution in the flow field. Sequences of full-order results are created for pressure 
contours to evaluated how the frequency affects pressure. Figs. 4.4 and 4.5 show the non¬ 
linear effect of frequency for three cases (0.1 in first column, 0.5 in second column, and 
1.0 in third column). The figures in each row show identical panel deflection amplitudes 
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Figure 4.3 Minimum (full-system) 

for difference frequencies, and each row is a different increment in time. As the panel is 
deflected in time the fluid is displaced in a wave that propagates away from the panel. 
At low oscillation frequencies (below 0.4), the problem is quasi-steady and the flow over 
the panel dominates the fluid structure. The pressure contours are nearly symmetric and 
there is little interaction upstream of the panel. As the oscillation frequency is increased, 
the velocity of the panel approaches that of the freestream fluid and pressure waves are 
created that interact upstream with the flow field to create non-linearities. The pressure 
waves are dissipated upstream by the inlet conditions and are passed downstream and out 
of the flow field. As the frequency increases further, the panel oscillates quickly enough to 
produce a second pressure wave before the first pressure wave can pass beyond the panel. 
This results in an amplification of the pressure waves as they combine upstream. 

4.2 POD/ROM Results 

Several experiments were completed using the oscillating panel model problem. First, 
comparisons were made between the accomplished accuracy of static grids as compared to 
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that of deforming grids. Next, POD/ROMs of various amplitudes and frequencies were 
created and tested to determine their accuracy. Then, these POD/ROMs were combined 
by merging their snapshot matrices to determine if a single model could be created. Fi¬ 
nally, comparisons were made of the accuracy of deforming grid POD/ROMs when run at 
amplitudes differing from their training amplitudes. 

4-2.1 Individual POD/ROM. To evaluate the accuracy of POD/ROMs with 
deforming grids, POD/ROMs were constructed for the static, POD/ROM/RG, and de¬ 
forming grid, POD/ROM/DG, cases. A single POD/ROM/RG was taken at a maxi¬ 
mum amplitudes (A) of 0.1 and an oscillation frequency, /, of 0.1. The low frequency 
resulted in a quasi-steady problem, necessary for the TBC to be accurate. Individual 
POD/ROM/DGs were constructed for deforming grids computed at maximum amplitudes 
{A = 0.05,0.1.0.15). At the various amplitudes, nine oscillation frequencies were used 
ranging from 0.1 to 1.0, creating 36 different POD/ROM/DGs. The POD/ROM/DG were 
trained for only one oscillation, regardless of frequency of oscillation. During the train¬ 
ing period, 20 snapshots were taken. The POD/ROM/DGs were then compared to their 
respective full-order solutions for fully developed flow (that occurred after approximately 
ten oscillations). 

Static Grid: For the quasi-steady frequency of 0.1 and amplitude of 0.1, the 
POD/ROM/RG was able to reproduce an accurate solution. The primary mode of the 
static grid case was a scalar fraction of the flow at the largest deflection and accounted 
for 96% of the total modal contribution (Fig. 4.6). The modes of the POD/ROM/RG 
show structures similar to wave shapes that propagate upstream (Fig. 4.6 (h), (i), and 
(1)). The remaining modes decreased in energy rapidly and were negligible zero by the 9th 
mode (10^^) (Fig. 4.7). The slope of the modal contribution was steep and indicated that 
relatively few modes were needed for an accurate solution. The POD/ROM/RG were able 
to reproduce the fluid flow at the fully developed condition (approximately 10 oscillations) 
very accurately with only 5 modes. This resulted in a reduction from 2400 to 5 dof. The 
flow field was reproduced with all fluid structures (Figs. 4.8 and 4.9). The results of the 
pressure contours for the POD/ROM/RG compared well to that of the full-order solution 
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(a) Mode 1 


(b) Mode 2 


(c) Mode 3 



(j) Mode 10 (k) Mode 11 (1) Mode 12 


Figure 4.6 POD/ROM/RG Modes 






Figure 4.7 POD/ROM/RG Modal Contribution 

(POD/ROM/RG in the first column and Full-order in second column) for the same fre¬ 
quency/amplitude combination. The contours of pressure were nearly identical, with only 
small variations in the far-held contours. 

Deforming Grid: The individual POD/ROM/DGs were able to reproduce the 
deforming grid how held for the different amplitude/frequency combinations using just 
15 modes to an accuracy of 3%, based on density. This equated to an order reduction 
from 2400 to 15 dof. As the problem became more non-linear in either amplitude or 
frequency, more modes were necessary for an accurate solution, as shown by the reduction 
in accuracy for the higher frequency and higher amplitude cases at 15 modes (Fig. 4.10). 
The accuracy of the more non-linear regions could be improved by using additional modes. 
The how held was reproduced with all huid structures very accurately. The results of 
the pressure contours for the POD/ROM/DG compared to that of the full-order solution 
for various frequencies showed good accuracy with only slight deviations in the higher 
frequency POD/ROM/DG. (Figs. 4.11 to 4.16 for results of A = 0.1 and / = 0.1,0.5,1.0). 
POD/ROM/RG is shown in the hrst column and the full-order solution is shown in the 


4-10 




(a) POD/ROM, t = 5 


(b) Full System, t = 5 



(c) POD/ROM, t = 6 


(d) Full System, t = 6 



(e) POD/ROM, t = 7 


(f) Full System, t = 7 



(g) POD/ROM, t = 8 


(h) Full System, t = 8 


Figure 4.9 POD/ROM/RG vs. Full-Order Solution (part 2) (/ = 0.1, A = 0.1) 


□ Full-System (A = 0.05) 



Figure 4.10 Minimum Cjq (Individual POD/ROM/DG, 15 modes) 

second column for the identical frequency/amplitude combinations. In examining the 
0.1 amplitude cases, the modes of the POD/ROM/DG were significantly different from 
the POD/ROM/RG (Figs. 4.17 to 4.19 for results oi A = 0.1 and / = 0.1,0.5,1.0). The 
/ = 0.1 POD/ROM/DG primary mode was similar to that of the POD/ROM/RG, but the 
remaining modes differed in both near and far-held structures. The / = 0.5 and / = 1.0 
POD/ROM/DGs had modes that were unique, with no corresponding similarities to the 
POD/ROM/RG modes. The modal contribution for the individual POD/ROM/DGs 
showed the increasing non-linearity caused by the increasing frequency (Fig. 4.20). In the 
more linear cases (lower frequency), the lower order modes decreased rapidly in energy. 
In the more non-linear cases (higher frequency), the lower order modes did not decrease 
in energy as quickly. The primary mode of the / = 0.1 POD/ROM/DG contained 81% 
of the total energy as compared to 70% for the / = 0.5 POD/ROM/DG and 69% for the 
/ = 1.0 POD/ROM/DG. 

Static vs. Deforming Grid Comparison: The accuracy of both the static grid 
and deforming grid POD/ROMs was dependent on the number of modes retained, however, 
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(a) POD/ROM, t = 1 


(b) Full System, t = 1 



(c) POD/ROM, t = 2 


(d) Full System, t = 2 




(e) POD/ROM, t = 3 


(f) Full System, t = 3 



(g) POD/ROM, t = 4 (h) Full System, t = 4 

Figure 4.11 POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 1) (/ = 0.1, 
A = 0.1) 
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(c) POD/ROM, t = 6 


(d) Full System, t = 6 




(e) POD/ROM, t = 7 


(f) Full System, t = 7 



(g) POD/ROM, t = 8 


(h) Full System, t = 8 


Figure 4.12 POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 2) (/ = 0 
A = 0.1) 
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(c) POD/ROM, t = 0.4 


(d) Full System, t = 0.4 




(e) POD/ROM, t = 0.6 


(f) Full System, t = 0.6 




X 


(g) POD/ROM, t = 0.8 (h) Full System, t = 0.8 

Figure 4.13 POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 1) (/ 
A = 0.1) 
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(c) POD/ROM, t = 1.2 


(d) Full System, t = 1.2 



(e) POD/ROM, t = 1.4 


(f) Full System, t = 1.4 



(g) POD/ROM, t = 1.6 (h) Full System, t = 1.6 

Figure 4.14 POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 2) (/ 
A = 0.1) 
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(c) POD/ROM, t = 0.2 


(d) Full System, t = 0.2 




(e) POD/ROM, t = 0.3 


(f) Full System, t = 0.3 




X 


(g) POD/ROM, t = 0.4 (h) Full System, t = 0.4 

Figure 4.15 POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 1) (/ 
A = 0.1) 
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(g) POD/ROM, t = 0.8 


(h) Full System, t = 0.8 


Figure 4.16 


POD/ROM/DG vs. Full-Order Solution Pressure Contour (part 2) (/ 
A = 0.1) 
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(a) Mode 1 


(b) Mode 2 


(c) Mode 3 


(d) Mode 4 


(g) Mode 7 


(j) Mode 10 

Figure 4.19 



- 


(e) Mode 5 


(f) Mode 6 


- 


(h) Mode 8 


(i) Mode 9 


- 




.. , 


(k) Mode 11 (1) Mode 12 

POD/ROM/DG Modes (/ = 1.0, A = 0.1) 




Figure 4.20 POD/ROM/DG Modal Contribution 

the static grid was more accurate with fewer modes (Fig. 4.21). The additional modes 
were due to the relative motion of grid points in the deforming grid. In the deforming 
grid case, at maximum amplitude grid points were compressed to fit between the deflected 
panel and the far field. This brought more grid points relatively closer to the panel surface 
and the higher dynamics of the flow field. In the static grid, the grid points were at a 
fixed distance from the panel surface (Fig. 4.2). In the transpiration boundary condition 
case, the far field grid points were 2.25 length units from the panel surface at all times. In 
the deforming grid case, at maximum amplitude, the far field grid points were 2.15 length 
units from the panel surface. This relatively closer distance resulted in a wider range of 
fluid dynamics seen by the far held grid points in the deforming grid case. This effect was 
similarly repeated in the grid points throughout the deforming region. The result was an 
increase in modes necessary to model the full-order solution. 

Deforming Domain Evaluation: In model problems in which a static grid cannot 
be used, another method was needed to isolate the effect of deformation on the POD/ROM 
accuracy. By reducing the number of grid points that are allowed to deform, the accuracy 
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Figure 4.21 Maximum Density Error (A = 0.1, 0.1 hz) 


Table 4 

1.1 Effect on Error due to Grid Density 


# Moving Pts 

% Max Error (2 modes) 

% Max Error (10 modes) 

100% Deforming Domain 

600 

7.76 

3.52 

15% Deforming Domain 

100 

5.64 

3.37 

Static Domain 

0 

1.99 

1.67 


trend could be observed. Two domains were evaluated: 100% and 15%. The 0.1 amplitude 
and 0.1 frequency case was used with 2 modes and 10 modes. (Tbl. 4.1). Maximum 
density error was evaluated. 

As the number of modes increased, the difference between the two domains became 
insignificant (> 0.5%). However, with fewer modes, the domain did have a small effect on 
accuracy (a reduction of 2.0% error from the 100% deforming domain to the 15% deforming 
domain. When compared to the static grid case, the trend is similar. 


4-2.2 Blended POD/ROM. To evaluate the robustness of POD/ROM/DGs, 
snapshots were combined for a variety of frequencies to determine if a single POD/ROM/DG 
could be developed for all frequencies tested. Snapshots from an amplitude of 0.1 and 
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A Full-System (A = 0.1) 
o Individual POD (A = 0.1) 

[> Blended POD (Trained at A = 0.1, Run at A = 0.1) 



Figure 4.22 Minimum (Blended POD, trained at A = 0.1, run at A = 0.1, 15 modes) 

frequencies of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 1.0 were combined to create a 
POD/ROM/BIO, in a technique called blending. 20 snapshots were taken for each fre¬ 
quency through a single oscillation, resulting in 180 snapshots for the blending. The 
POD/ROM/BIO was then run at the same frequencies and compared to the respective 
full-order at fully developed flow. The POD/ROM/BIO was stable at all frequencies (Fig. 
4.22) but less accurate than the individual POD/ROM/DG, especially for the higher fre¬ 
quencies. The maximum density error for all frequencies was within 7%. The modes of the 
POD/ROM/BIO were again unique. The higher energy modes contained wave structures 
that appeared to be combinations of the wave forms seen in the individual POD/ROM/DG 
(Fig. 4.23). The primary mode contained 63% of the total eigenvalue energy. The remain¬ 
ing modes decreased in a fashion similar to the higher frequency individual POD/ROM/DG 
(Fig. 4.24). However, the slope of the POD/ROM/DBIO modal contribution was shallower 
than that of the individual POD/ROM/DGs, indicating more modes were necessary for 
an accurate solution. 
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(a) Mode 1 


(b) Mode 2 


(c) Mode 3 



(j) Mode 10 (k) Mode 11 (1) Mode 12 


Figure 4.23 POD/ROM/BIO (Blended, A = 1.0) 




Figure 4.24 POD/ROM/BIO Modal Contribution 

4-2.3 Amplitude Variation. To evaluate changes in deformation on the accuracy 
of deforming grids, the blended POD/ROM trained at an amplitude of 0.1 was run at 
amplitudes of 0.05 and 0.15. This provided a ±50% variation in the grid deformation. 
The blended POD/ROM run at the lower amplitude of 0.05 was only slightly less accurate 
than that run at the amplitude at which it was trained (Fig. 4.25). The maximum density 
error was 5% for all frequencies. The blended POD/ROM run at the higher amplitude 
of 0.15 was much less accurate than that run at the amplitude for which it was trained 
(Fig. 4.26). The maximum density error for all frequencies was ±12%. These amplitude 
variation errors were due to the modes being applied at inappropriate physical locations. 
To created a single model, a blended POD/ROM using the snapshots from 0.05, 0.1, and 
0.15 with all frequencies, was developed. This single model was unstable at all amplitudes 
and frequencies due to the fact that the modes combined across amplitudes were being 
applied at inappropriate physical locations. This highlights the issue with deforming grid 
POD/ROM; snapshots of differing deforming grids cannot be easily combined to produce 
a blended POD/ROM. As a result some new technique was required. 
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Figure 4.25 


Figure 4.26 



Frequency 


Minimum Cn (Blended POD, trained at A = 0.1, run at A = 0.05, 15 
modes) 



Frequency 


Minimum Cn (Blended POD, trained at A = 0.1, run at A = 0.15, 15 
modes) 
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4.3 MulU-POD 

To create a more robust model, a technique that used multiple POD/ROM/B was 
developed (Multi-POD). In this technique, several independent POD/ROM/DG were de¬ 
veloped at various sets of grid deformations. Blended POD/ROM were then created by 
combining snapshots for the POD/ROM/DG. POD/ROM/B tend to be more effective than 
individual POD/ROM/DG because they are more robust to variations in input parame¬ 
ters. The accuracy of the Multi-POD is dependant on the accuracy of the POD/ROM/B. 
Blending was not effective over variations in grid deformations due to the changes in the 
grid point locations caused by the grid deformation, as seen in the pervious section. Using 
the deforming grid metric of the L 2 error norm, several POD/ROM/B were used for a 
single model that had a variety of deformations. As the deformation progressed, the L 2 er¬ 
ror norm deforming grid metric was automatically checked. When a user-defined level was 
exceeded, the algorithm automatically switched to a new POD/ROM/B, more appropriate 
to the current deformations. 

As an example, the amplitude variation errors in the POD/ROM/B 10, created in 
the previous section, when it was run at different amplitudes could be correlated to the 
differences in the training grids to the run grids. The L 2 norm of the grid error was well 
correlated to the POD/ROM/BIO error (Fig. 4.27). As the run grids deviated from the 
training grids, the overall POD/ROM/BIO accuracy decreased, as shown by the L 2 norm 
leading the density error. The L 2 norm was scaled by 100 to allow a visual comparison. 
For this model problem, L 2 norms of the grid error of around 0.065 produced maximum 
density errors of 8.5%. As a comparison, the transpiration boundary condition case would 
have a grid error of zero at all times, providing no chance to correlate the errors. The L 2 
norm of the grid error was provided to the user as a flag on the accuracy of the POD/ROM 
for the deforming grid. 

To test the Multi-POD technique, two blended POD/ROMs, trained at amplitudes 
of 0.05 and 0.15 (POD/ROM/B5 and POD/ROM/B15), were used to model amplitude 
deffections from 0.05 to 0.15. In this example, the problem was run at an amplitude of 
0.05 for five oscillations, then 0.15 for five oscillations. Starting with POD/ROM/B5, 
the model proceeded up to an amplitude of approximately 0.1. At this point, the L 2 error 
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Density Error (%) 

L2 Norm of Grid Error (x100) 



Figure 4.27 Correlated Grid Error and Density Error, after 8 oscillations (Blended POD, 
trained at A = 0.1, run at A = 0.15, 15 modes) 

norm exceeded the preset level and POD/ROM/B15 took over. When compared to a single 
blended POD/ROM (trained at an amplitude of 0.1), the Multi-POD was more accurate 
(Fig. 4.28). The density error of the POD/ROM/BIO exceeded 8% while the Multi-POD 
remained below 5%. The L 2 error norm clearly identified the point at which a change 
in POD/ROM was necessary as shown by the sudden increase in the L 2 error norm plot. 
This technique required more snapshots than the single blended POD/ROM (double in 
this example). However, it does provide for a more robust algorithm that can be applied 
in a design environment. With the multi-POD technique, the entire design space for the 
example is well modeled with a single algorithm containing two POD/ROMs. 

4-4 Conclusions 

POD/ROM was shown to be sensitive to grid deformation. When compared to a 
rigidly attached grid, the deforming grid required more modes for an accurate solution. 
The relative motion of grid points due to the deformation created modes that did not exist 
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L2 Grid Error Norm (Blended POD) 



Figure 4.28 Multi-POD (Two POD; trained at A = 0.1 and A = 0.15, 15 modes) 

in the fluid dynamics. The index-based computation of POD/ROM could not account for 
such changes in grid dynamics. 

POD/ROM was shown to accurately reproduce flow solutions with deforming grids 
for a problem with coupled fluid and grid dynamics. Errors in fluid variables were less 
than 5% of the full-order system, and thus represent a reasonable level of accuracy for 
design. POD/ROM accurately reproduced the flow field when applied to a deforming grid 
model problem when the run grid was identical to the training grid. The maximum density 
error was 3% for only 15 modes. POD/ROMs of deforming grids required more modes for 
accuracy than those of static grids, for similar model problems. The additional modes 
were attributed to the relative motion of the grid points with respect to the boundaries 
and each other. However, the overall number of modes required for the deforming grid was 
still relatively small, reducing the order of the problem from 2400 to 15 dof. 

The magnitude of the deforming domain provided a technique for evaluating the 
accuracy of POD/ROM on deforming grids. As the domain was reduced, the accuracy 
increased for a similar number of modes. The trend became constant as the domain reduced. 
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indicating the number of modes that were associated primarily with grid deformation. This 
constant was similar to that of the static grid case. 

When the modes were applied at grids that differed from the training grids, POD/ROM 
accuracy degraded. A POD/ROM trained at 0.1 amplitude and run at 0.15 amplitude had 
errors of ±12% in maximum density. A metric was developed to track the relative differ¬ 
ence in training grids and run grids. A pseudo-inverse technique was implemented from a 
grid modal matrix based on the training grids and an error term created. The ±2 norm of 
the error in deformation correlated well with the resulting POD/ROM error. The metric 
was used interactively to evaluate the accuracy of the deforming grid POD/ROM. A new 
multi-POD technique was developed to use the grid metric to determine when to switch 
between POD/ROM of different deforming grids. The complete algorithm used the most 
appropriate POD/ROM for a variety of grid deformations to provide the most accurate 
solution. 
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V. Pitching and Plunging Airfoil 

To examine the utility of the Multi-POD technique with a deforming grid, a pitching and 
plunging airfoil with a loosely coupled structural model governing the pitch and plunge 
degrees of freedom was analyzed. The structural model allowed the aerodynamics to alter 
the grid deformation by movement of the structure through imposed air loads. This free 
pitching and plunging required that the Multi-POD technique be capable of selecting the 
best available POD/ROM for cases where the pitch and plunge were not well represented 
in the training. The training of the POD/ROM was performed without the structural 
model, using cases for which pitch and plunge was forced. 


5 .1 Application 

The application of full-potential flow over a 2-D pitching and plunging airfoil was 
used to examine the robustness of the Multi-POD technique (Fig. 5.1). Two cases were 
examined: forced oscillation and free pitch and plunge. For both cases, a NACA 0012 
airfoil was used. The airfoil was pitched about the quarter chord point. 

The density and velocity were non-dimensionalized using the far-field density, 
and velocity, Uoo- Pressure was non-dimensionalized using the far-field dynamic pressure, 
Poo (subsequent references to variables assume non-dimensional forms). For both cases, the 
freestream Mach number was set to 0.5 to prevent the flow over the airfoil from becoming 
supersonic. 


5.1.1 Forced Oscillation Case. In the Forced Oscillation case, both plunge, h {t), 
and pitch, a (t), were varied in a sinusoidal fashion; 


a (t) 
h{t) 


[cXq CX 

max) 


sin 



sin 



(5.1) 

(5.2) 


ao was the initial conditions, amax and hmax were the maximum pitch angle and plunge, / 
was the non-dimensional frequency of oscillation, t was the time, and tmax the maximum 
time for the test. The pitching and plunging was synchronous, with the pitch angle zero 
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Figure 5.1 Pitching and Plunging Airfoil Model Problem 


'■ ! 



Figure 5.2 Pitching and Plunging 2-D System 

at minimum and maximum deflection. The grid was deformed by determining the pitch 
and plunge of the airfoil and then analytically solving for the grid-points from the airfoil 
surface to the fixed far-held. For the forced case, the pitch and plunge was explicitly known 
and not dependant on the huid dynamics. 

5.1.2 Free Pitch and Plunge Case. For the free pitch and plunge case, the 
structural model was a simple, 2-D, spring model that was loosely coupled to the huid 
dynamics model (Fig. 5.2). The huid dynamics model generated the coefficients of lift 
and moment. The structural model then used the huid dynamics loads to determine 
incremental motion of the airfoil. This model allowed an evaluation of the POD/ROM 
with the structural dynamics and grid motion fully linked. The airfoil was free to rotate in 
the x—y plane and free to translate up or down, about the quarter chord point. The motion 
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of the airfoil was limited by spring effects in torsional and translational displacement. The 
equations governing the motion of the system are, 


h H—T + C2/1 — - Cl 

2 

XqJi H—H-+ Cl (a — ao) = - 

2 U HgTT 


where, ci = ^ *^2 = (|) h = 5 ^, and Xa is the distance from the 

center of mass to the pitch point, scaled by the semi-chord. The pitch and plunge natural 
frequencies are defined as, zua = ^^ and zuh = The linear and torsional springs 

are modeled with stiffness coefficients, and Ka, and damping coefficients, and Da. 
The static pre-twist, ao) defines the unloaded position of the torsional spring. The radius 


of gyration is, = y —. The three remaining parameters are the mass ratio, , 

and the damping ratios, Cq, = 2 Ji°‘k ~ 2 jmK ' These equations can be rewritten 

as a system of first-order differential equations by setting yi = h, y 2 = h, 2/3 = a, and 
2/4 = a. Following the approach of Beran and Morton (42), with det = ^ (r^ — , 


S = /(S,a,C'^..,h) = M-iQ-M-iKS. 



det 0 0 0 

1 0 0 -^Xa 

0 0 det 0 

0 —Xa 0 1 


(5.5) 


(5.6) 


(5.7) 
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0-10 0 

4 2CftC2 0 0 

0 0 0 -1 
0 0 Cl Ca(i)c« 

The above equations may then be then solved using a four-stage Runge-Kutta method. 
In the structural model, the parameters used were selected to closely match the forced 
oscillation case in terms of amplitude and frequency (Tbl. 5.1). 

The grid was deformed based on the pitch and plunge generated from the above 
equations and the air loads calculated from the fluid dynamics model. This implicitly 
linked the fluid dynamics and grid motion. The structural solver was validated using 
computational data (see Appendix D). 

5.1.3 Computational Grid. For both cases, a C-grid was used with a wake cut 
along the trailing edge that extended to the outflow boundary (Fig. 5.3). The trailing 
edge of the airfoil was closed. The chord length, c, was normalized to a value of 1.0. The 
far field was placed at 10 chord lengths from the airfoil. 251 x 50 grid points were used 
for the model problem, with 150 grid points spaced on the panel surface. The spacing of 
the grid at the airfoil surface was 0.05 of the chord length. The grid was smoothed with 
an elliptical solver that used a tension spline to prevent grid line crossing (Fig. 5.4). 

5 . 1.4 Domain Filtering. In POD/ROM, the most accurate solution is typically 
generated by the most independent modes. However, with a large region of nearly constant 
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(a) Computational Domain (b) Grid Region Around Airfoil 

Figure 5.3 C-Grid 

flow, large numbers of modes can generate numerical errors. An approach to dealing with 
this numerical error is to sub-divide the computational domain. For the regions of highest 
interest (near field), more modes were used. For areas of lower interest (intermediate field), 
fewer modes were needed. In areas of little change from the freestream (far field), only a 
few modes were necessary (Fig. 5.5). 

To implement this, a domain-based filtering of the POD modal matrix, d>, was per¬ 
formed. For the grid-point locations to be filtered, the value of $ corresponding to the 
appropriate modes, was set to zero. This truncated the number of modes being applied to 
various regions of the computational domain. In the near field, the number of modes var¬ 
ied depending on the complexity of the fluid flow (see POD/ROM Results for numbers of 
modes). For both the intermediate and far fields, it was determined through trial and error 
the number of modes necessary for an accurate solution. If there was too large a variation 
in number of modes from the near field to the intermediate field, errors would develop 
along the internal boundary. In the intermediate field, only 10 modes were necessary. For 
the far field, only 2 modes were necessary. This domain filtering technique differs from the 
Domain Decomposition technique of Lucia et al. (43), in that only one POD/ROM was 
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Figure 5.4 Example of Grid Deformation (/imax 


0 .1, Umax = 2.0 deg) 





































Far Field 



Figure 5.5 Grid Domains 

generated for the entire domain. However, the technique of domain filtering required more 
computations during the actual run time. 

Other forms of filtering the domain include weighted filtering. In weighted filtering, 
the modes are multiplied by weighting factors that decrease with radial distance from the 
airfoil. In effect, the domain filtering technique is a binary weighted filtering. Weighted 
filtering could be explored but was not necessary for this model problem. 

5.1.5 Full-Order Results for Forced Oscillation. A wide range of input parameters 
was evaluated. Maximum pitch was varied from 0 to 2.0 degrees in 0.5 degree increments. 
Maximum plunge was varied from 0 to ±0.2 of chord length, in 0.05 increments. Oscillation 
frequency, /, was varied from 0.002 to 0.01, in 0.002 increments. This provided 125 different 
combinations of input parameters. The model problem showed non-linear behavior in both 
maximum plunge and maximum pitch. The maximum lift coefficient, Ci, determined after 
the flow was fully developed, was used as a way to evaluate the non-linear nature of 
the model problem. Examining the case with maximum plunge of 0.2 and varying pitch 
magnitude and frequency, the model problem can be seen to be non-linear in pitch for 
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the lower frequencies (Fig. 5.6). As the oscillation frequency increased, the plunge rate 
became more significant, overshadowing the pitch. Examining the case with maximum 
pitch of 2.0 degrees and varying plunge and frequency, the model problem can be seen to 
be non-linear in plunge as well. The rate of plunge causes the apparent angle of attack to 
be increased with respect to the freestream flow (Fig. 5.7). 

5.2 POD/ROM Results 

POD/ROMs were created from the forced oscillation cases and then applied to both 
the forced and free dynamics. Individual and blended POD/ROMs were compared to the 
forced oscillation case. Then, blended POD/ROMs and Multi-POD were compared to the 
free pitch and plunge case. 

5.2.1 Individual POD/ROM. Individual POD/ROMs were constructed for the 
full-order solutions computed at pitch angles, amax = 0.0 — ±2.0 deg in 0.5 deg increments 
and plunge depths /imax = 0.0 —±0.2 in increments of 0.05 for oscillation frequencies of / = 
0.02 — 0.1 in increments of 0.02. These parameters gave a varied design space with sufficient 
non-linear variation in behavior to rigorously test POD/ROM with deforming grids. A total 
of 125 separate POD/ROMs were developed. The individual POD/ROMs were trained for 
two oscillations and 50 snapshots were taken during the training period. The POD/ROMs 
were then compared to their respective full-order solutions for fully developed flow (that 
occurred after approximately 10 oscillations). 

The individual POD/ROMs were able to reproduce the flow field using just 20 modes 
to an accuracy of 5%, based on density throughout the flow-field. This equated to an order 
reduction from 12550 to 20 dof. The POD/ROMs were able to capture the non-linearities 
in both pitch and plunge with excellent accuracy (Figs. 5.8 and 5.9). The lift coefficient was 
reproduced to an accuracy of 8%. The primary modes of all of the individual POD/ROMs 
were very similar in magnitude, however, the modal shapes were very different. For low 
frequency cases (/ = 0.02 and 0.04) the modes were simpler with little data in the near 
and far-flelds. For the higher frequency cases, the modes showed more information in the 
far-fleld. In the case of varying pitch, frequency of plunge determines the mode shapes. For 
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Figure 5.6 Full-Order (Fixed Plunge, Varying Pitch) 
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Figure 5.7 Full-Order (Varying Plunge, Fixed Pitch) 








Figure 5.8 Individual POD/ROM (Fixed Plunge, Varing Pitch, 20 modes) 




lower frequencies the modes were similar. However, as the frequency increases the modes 
changed significantly (Figs. 5.10 and 5.11). A similar result was seen in the frequency of 
pitch. The pitch angle had a less significant impact on the modal shapes. (See Appendix 
E. for additional modal shape comparisons). The first three modes for all cases had 
similar magnitude modes (Fig 5.12). The modal contribution for the remaining modes 
depended on the complexity of the modes. If the modes were less complex, such as in a 
lower frequency pitch and plunge case, the majority of the modal energy was contained in 
the primary mode. However, if the modes were more complex, the modal energy is spread 
out over more modes. 

5.2.2 Blended POD/ROM. Blended POD/ROMs were constructed from the 
individual POD/ROM for plunge depths of hmax = 0.0, POD/ROM/BO, /imax = 0.1, 
POD/ROM/BIO, and /imax = 0.2, POD/ROM/B20. The snapshots for the various os¬ 
cillation frequencies and pitch angles were combined for each of the three plunge depths. 
The blended POD/ROMs were then run at the same frequencies and compared to the re¬ 
spective full-order solution at fully developed flow. The maximum flow-field density error 
of the blended POD/ROMs over their range of oscillation frequency and pitch angle, was 
within 6% and the lift coefficient accuracy was within 10% using 20 modes (Fig 5.13). 

If a blended POD/ROM was run on a grid that did not match it’s training grid, the 
solution was very poor. The results of POD/ROM/BO when run at /imax = 0.1, / = 0.06, 
and a = 1.0 deg were very unstable (Fig 5.14). POD/ROM/BO was unable to accurately 
reproduce the non-linear behavior of the full-order system and became more unstable as 
the frequency was increased. The maximum Ci was damped for all frequencies and shifted 
dramatically as a result of the modes being applied at inappropriate physical locations. 

The accuracy of the blended POD/ROM depended on the number of modes used. For 
the more quasi-steady cases, fewer modes could be used to obtain a reasonable accuracy. 
In the more non-linear cases, more modes were necessary. POD/ROM/BO, required 10 
modes for an accurate solution. POD/ROM/BIO required 15 modes and POD/ROM/B20 
needed 25 modes (Fig 5.15). The blended POD/ROM showed significant structures in the 
far-held where the huid dynamics was in fact quite small. This was due to blended over 
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Figure 5.12 Individual POD/ROM Modal Contributions 
Table 5.2 Computation Times vs. Number of Modes 



100 Iterations (sec) 

full-order 

100 

5 Modes 

53 

10 Modes 

71 

15 Modes 

94 

20 Modes 

128 

25 Modes 

145 

30 Modes 

223 


the higher frequency cases that had modes with data in the far-held. While the individual 
POD/ROM were able to cancel out the modes in the far-held, the blended POD/ROM were 
not. These blended far-held structures resulted in increasing numerical error in the higher 
modes. Thus, a practical limit was seen in the number of modes that could be used in a 
POD/ROM. This was seen when increasing the number of modes actually increased the 
error. The fewer the modes used in the blended POD/ROMs, the faster the computation 
of the reduced order model. Tbl. 5.2 shows the computation time for 100 iterations with 
5 sub-iterations, on a Pentium III, 933 Mhz, with 512 MB of memory. Therefore, it was 
more computationally efficient to use blended POD/ROM with the appropriate number of 
modes. 


5-15 




Maximum Cl Maximum Ci 


^ - Full-Order 




(a) /imax = 20% 


(b) /ln,ax = 10% 



(c) /imax = 0% 

Figure 5.13 Blended POD/ROM (Fixed Plunge, Varying Pitch, 25 Modes) 
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Full-Order 



Figure 5.14 Blended POD/ROM (Trained at /i = 0%, Run at h = 10%, / = 0.06, 
a = 1.0 deg, 25 modes) 



Figure 5.15 Blended POD/ROM Accuracy vs. Number of Modes 
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The modes of the blended POD/ROMs were more complex than that of the individual 
POD/ROMs, with large structures in the first eight modes (see Figs. 5.16 to 5.18). The 
modes showed more structure in the far field. As a result of the large structures in the 
far-field, the POD/ROM/B were more sensitive than the individual POD/ROM/DG to 
numerical errors caused by using too many modes. Thus, domain filtering was vital for an 
accurate solution using blended POD/ROM. The energy of the modes was concentrated 
in all cases in the first four modes (see Fig. 5.19). For POD/ROM/B20, the remaining 
modes seemed to drop off quickly in energy content, however, the modes were still vital 
to an accurate solution. This is due to the fact that more energy is contained in the first 
four modes than in any other case. For POD/ROM/B 10, the modes dropped off in energy 
more slowly, more like in the case of the individual POD/ROMs. 

5.2.3 Multi-POD. To determine the utility of Multi-POD for a practical problem, 
the blended POD/ROM created using forced pitch and plunge were applied to the free pitch 
and plunge cases. The three blended POD/ROMs were assigned thresholds for the L 2 of 
the deforming grid metric. 10 Modes were used for POD/ROM/BO, 15 for POD/ROM/BIO 
and 20 for POD/ROM/B20. 

For the damped case, the model was run with u = 2.0 and the Multi-POD compared 
to the full-order solution. The Multi-POD switched to POD/ROM/BO very quickly as 
the airfoil came down from a = 1.0 deg and then held there until it was fully damped 
(see Fig. 5.20). The Multi-POD was able to damp the solution but did not match the 
time accurate response in every detail. The phase of the motion of the pitch and plunge 
matched well, but the amplitudes of both pitch and plunge were larger than that of the 
full order solution. The Multi-POD was exact at steady state (fully damped). To test 
how well a single POD/ROM would have performed, the same test was executed with the 
POD/ROM/BIO. The solution was less accurate during the unsteady portion of the time 
dependent solution, but exact at steady state (fully damped). 

The unstable oscillation case was more challenging. The model was run with u = 5.0 
and the Multi-POD compared to the full-order solution. The Multi-POD started with 
the POD/ROM/BO and switched as the amplitude varied slowly to POD/ROM/BIO. The 
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Figure 5.16 Blended POD/ROM Modes (/imax 
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Figure 5.21 Multi-POD, Unstable Case {u = 5.0) 

Multi-POD varied between the two POD/ROMs until it finally switched to POD/ROM/20 
as the solution became unstable. The solution was damped, but the non-linear behavior 
of the full system was accurately modeled (Fig 5.21). The changes between POD/ROMs 
were smooth and did not materially affect the solution. The Multi-POD accuracy was 
within 5% for pitch and 10% for plunge. To test how well a single POD/ROM would have 
performed, the same test was executed with the POD/ROM/BIO solely. The solution was 
heavily damped and did not match the full-order solution in a time accurate sense (see 
Fig 5.22 for variations). The errors were significant in pitch (over 100%), and could not 
capture the dynamics of the full-order solution. 

5.3 Conclusions 

POD/ROM was shown to be effective for problems with complex, coupled structural 
and fluid dynamics models. Errors in fluid variables were less than 6% of the full-order 
system for this model problem, and thus represent a reasonable level of accuracy for design. 
In the dynamic case of a pitching and plunging airfoil, POD/ROM was able to reproduce 
the full-order solution to within 6%. Individual POD/ROMs required only 20 modes for 
an accurate solution, producing a three order of magnitude reduction in dof (12250 to 20). 
Blended POD/ROM combined a variety of frequency and pitch angles into a single model 
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Figure 5.22 Blended POD (Trained at ^max = 10%), Unstable Case {u = 5.0) 

that accurately described a single set of grid deformations. Blended POD/ROMs required 
slightly more modes for accurate solutions in the higher dynamic cases. 

The Multi-POD technique was shown to be effective for practical fluid dynamics 
problems by use of the model problem. The Multi-POD technique was able to reproduce 
the free deformation solution very accurately, based on POD/ROM that were developed 
from forced deformation cases. The Multi-POD was able to switch to the best available 
POD/ROMs for the deformation created by the loosely coupled structural model. 
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VI. Summary 

The thesis of this dissertation, POD/ROM can be used to simulate fluid flows accurately 
on deforming grids, has been demonstrated through analysis and application to example 
problems. In POD/ROM, snapshots are taken of the unsteady flow-field data at inter¬ 
vals of time during a process called training. These snapshots are used to generate fluid 
modes. The modes are then used to create a reduced-order model. If the snapshots repre¬ 
sent the fluid dynamics of the problem sufficiently well, the POD/ROM will be accurate. 
POD/ROM computations are conducted in computational (or index) space and do not di¬ 
rectly account for grid deformation. The order of the data in the index space is irrelevant 
to the process, but must be consistent. At each snapshot in time, the grid is assumed to be 
identical to all other grids. In a problem with a deforming grid, the physical location of any 
data point changes from one snapshot to the next, even though its location in index space 
does not. Therefore, the set of grid point locations captured by the snapshots represent 
the range of grid deformation for which the problem was trained. In the simplest terms, 
variations in grids cause variations in modes created by POD/ROM. If two deforming grids 
are sufficiently different, the modes from one POD/ROM cannot be used with the other 
grid. 

POD/ROM was shown to be sensitive to grid deformation. When compared to a 
rigidly attached grid for a nearly identical model problem, the deforming grid required more 
modes for an accurate solution. This shows conclusively that deforming grid POD/ROM 
differ from rigid grid POD/ROM due to the deformation. The additional modes necessary 
for accuracy for the POD/ROM/DG on a nearly identical problem were attributed to the 
relative motion of grid points due to the deformation. As the amount of deformation 
increased, in terms of the number of moving grid points, the number of modes necessary 
for an accurate solution increased. 

POD/ROM accurately reproduced the flow field when applied to a deforming grid 
model problem when the run grid was identical to the training grid. The overall number 
of modes required for the deforming grid was still relatively small, reducing the order of 
the problem by three orders of magnitude. When the modes were applied at run grids that 
differed from the training grids, the POD/ROM accuracy degraded. 
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A technique was developed that evaluated the relative difference in training grids 
and run grids. The L 2 norm of the error in deformation correlated well with the resulting 
POD/ROM error. The metric was used in a real-time evaluation of the accuracy of the 
deforming grid POD/ROM, in a technique called Multi-POD. The Multi-POD technique 
determined the current relative grid deformation and selected the best POD/ROMs from 
those available. This resulted in a more robust algorithm that could be used over a variety 
of deformations and input parameters. 

Finally, Multi-POD was applied to a problem where the grid deformation varied based 
on the fluid dynamics. In this case, separate POD/ROM were trained with forced grid 
deformation. The number of modes used in each POD/ROM was based on the complexity 
of the flow, with simpler flow fields requiring fewer modes. Then the model was run 
with free grid deformation. The Multi-POD technique was able to reproduce the free 
deformation solution very accurately, switching between POD/ROM as necessary. This 
demonstrated the robustness of the overall technique to a realistic model problem. 

6.1 Significant Advances 

Three significant advances were made in this work in answering the research questions 
proposed in the Introduction. 

6.1.1 Static vs. Deforming Grid POD/ROM Analysis. In analyzing the amount 
of error associated with POD/ROM/DG, a technique was developed to compare various 
amounts of grid deformation by varying the number of deforming grid points. The sensitiv¬ 
ity of POD/ROM accuracy to grid deformation is problem dependent. The degradation in 
accuracy due to grid deformation is based on the extent, relative magnitude, and scale of 
the deformation. By comparing the accuracy of the POD/ROM/DG created with varying 
numbers of deforming grid points for the same model problem and number of modes, a 
trend can be developed that shows the expected decrease in accuracy of the POD/ROM 
due to grid deformation. 
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6.1.2 Multi-POD. The Multi-POD approach is a unique technique that allows 
the best available POD/ROM to be used for a range of grid deformations, by switching 
between various POD/ROMs based on the deforming grid metric. A metric was developed 
to compare training and application grid deformation. This makes POD/ROM more robust 
over a wider range of grid deformation, and consequently over a wider set of applications. 
Without Multi-POD, POD/ROM/DG would be of limited utility. 

6.1.3 Application to Aeroelastic Problem. A unique application was executed 
using POD/ROM/DG. POD/ROM/DGs were trained using a forced grid deformation and 
then applied to free grid deformation, where the grid dynamics were fully coupled with 
the fluid dynamics. The training of a POD/ROM/DG in a forced case is more practical, 
and allows the user to define a domain of solutions of interest. Then, in applying the 
POD/ROM/DGs to a free case, a wider range of problems of interest can be evaluated. 

6.2 Future Work 

There are three directions for proposed future work. First, implement the Domain 
Decomposition technique of Lucia, et. al. (43). In most practical problems, the deforming 
grid is restricted to a relatively small region around the body of interest. The domain de¬ 
composition technique allows the creation of a large POD/ROM for the deforming domain 
and smaller POD/ROM for the static grid domains further out from the body. The static 
grid domains would require fewer modes than the deforming grid domains, due both to 
the lack of deformation and the less significant fluid dynamics in the far field. Also, only 
the deforming domain would necessarily require a Multi-POD application. This would 
significantly save computation time. 

Second, the full-order solver should be changed to a Navier-Stokes viscous solver. 
The full-potential solver is not robust enough to tackle complex model problems and the 
effect on POD/ROM modes due to viscosity in the deforming region should be explored. 
A new solver would also allow for an in-depth analysis of the computational savings that 
could be expected by using POD/ROM with deforming grid problems. The full-potential 
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solver uses a tri-diagonal matrix inversion that cannot be duplicated by POD/ROM and 
therefore any computational savings analysis is skewed in favor of the full-order solution. 

Finally, the grid dynamics could be incorporated into the POD analysis. The x 
and y grid point calculations could be treated as fluid variables and incorporated into the 
implicit Jacobian matrix, allowing POD to be applied to them. In addition, the structural 
model could be incorporated in the Jacobian matrix, allowing for a fully coupled aeroelastic 
model. These two changes would allow the POD/ROM to be more robust to changes in 
the grid deformation, although not eliminating the need for Multi-POD. 
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Appendix A. Potential Flow Around a Translating Cylinder 

To determine the impact of grid deformation on POD/ROM, an attempt was made to 
isolate grid deformation effects from a numerical solution of fluid flow. A problem with an 
analytical solution was selected: potential flow around a cylinder. This allowed the fluid 
dynamics to be de-coupled from the grid dynamics. The nature and complexity of the 
grid modes could then be analyzed to see how they might affect the numerical solution. 
Two grid cases were selected. The first case was a rigid grid fixed to the cylinder and 
translating with it. This case included no deformation effects. The second case had a 
fixed outer boundary and the cylinder translated within the boundary, deforming the grid. 
For each grid case, POD/ROM was applied to snapshots of data obtained by imposing 
the analytical solution on the grids at selected increments in time. With the analytical 
solution, difference in the POD/ROM of the two cases could be solely attributed to the 
grid deformation. 

A.l Model Problem 

Two grids were used in the model problem. In the first grid, the cylinder and grid 
rigidly translated together (Fig A.l). This eliminated any relative motion between the 
grid points. The cylinder oscillated through still fluid in the x direction and remained 
constant in the y direction. The center of the cylinder oscillated from —5.0 < x < 5.0 
(non-dimensional units). This was done by varying the velocity of the center of the cylinder 
using a cosine function. Time was divided into 100 units (0 < t < 100). In this grid, at 
every time increment, each grid point had the same relative distance to the cylinder. 

The second grid had an exterior boundary of fixed location, and the cylinder trans¬ 
lated within it. The grid deformation was accomplished in three ways (Fig. A.2). First, a 
simple translation of the cylinder. Second, a rigidly translating grid and outer boundary 
that rotated. Third, a translating of the cylinder with a stationary in translation but 
rotating outer boundary. The three sub-cases of the deforming grid were intended to sim¬ 
ulate the grid response of typical airfoil motions; plunging (translating), angle of attack 
(rotation), and combined pitching and plunging (translating/rotating). As in the first grid 


A-1 




(a) Rigidly Attached Grid (b) Deforming Grid 

Figure A.l Grid Cases for Translating Cylinder 

case, the cylinder was oscillated from —5.0 < x < 5.0 (non-dimensional units) by means 
of varying the velocity of the cylinder. When the outer boundary was rotated, the angle 
between the outer grid points and the corresponding cylinder points were evenly divided 
and the intervening grid points rotated incrementally. This resulted in a clocked grid. The 
cylinder was translated within the grid, causing deformation. This caused compression in 
the grid as well as rotation relative to the cylinder. The outer boundary is rotated from 
about —30° < a < 30°, 

a (t) = bcos {6t). (A.l) 

Time was divided into 100 units (0 < t < 100). The cylinder, with radius R = 1, was 
translated in the x direction at varying velocity, U (t), through a fluid otherwise at rest, 
which resulted in the perturbation stream function (44) 


T 

U{t) 


U (t) 
a sin {Ot) 


{x-xq {t)f + y‘^' 


(A.2) 

(A.3) 


where 4/ was the stream function and xq (t) the center of the cylinder at time t. 
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(a) Translating Grid Sub-Case 


(b) Rotating Grid Sub-Case (c) Translating/Rotating Grid 

Sub-Case 


Figure A.2 Grid Sub-Cases for Deformation 


To conduct an error analysis of POD/ROM with the deforming grid, snapshots were 
taken of the solution at various grid positions and cylinder velocities. This recorded data 
made up the snapshot matrix and was used to generate the modal matrix, 




(A.4) 


Next, a pseudo-inverse of $ was needed to compute T directly. This was accomplished using 
the Generalized Inverse (38), by multiplying both sides of Eqn. A.4 by the transformed 
matrix, This produced a square matrix, on the RHS of Eqn. A.4. By assuming 
that the approximation of 4/ is exact, was inverted and T determined. Einally, the 
error of the POD/ROM was determined by substituting T into the original approximation 
for 4/: 


errorq, = T — 4>T, (A.5) 

error^f = 4/ — 4> ^ 4>^4/. (A.6) 

The error resulted solely from the loss of information that occurred as the POD/ROM 
reduced the analytical solution to a smaller subspace. If $ were square and nonsingular 
(the number of modes equaling the number of grid points), the error in equation (A.5) 
would vanish and the POD/ROM would be exact. 


A-3 



Figure A.3 POD/ROM/RG vs. Analytical Solution (Maximum value of 'I/, lines are 

identical) 

A. 2 Results 

The POD/ROMs of the two grids were significantly different. The rigid grid case 
POD/ROM, POD/ROM/RG, was an exact solution and required only a single mode. 
The deforming grid POD/ROM, POD/ROM/DG, was not an exact solution and required 
modes for accuracy. The additional modes of the deforming grid were attributed solely to 
the deforming grid and considered grid modes (modes that represent the deformation of 
the grid and the subsequent mapping of the undeformed solution to the deformed space). 

A.2.1 Rigid Grid Case. The first mode of the POD/ROM/RG contained 100% 
of the total energy of the eigenvalues. The remaining modes were machine zero (10^^^). 
By plotting the contribution to the maximum value of the stream function for each mode, 
the relative contribution of each mode could be examined. The first mode and the true 
solution were identical lines (Fig. A.3). The single mode, when plotted against the rigid 
grid, was a scalar fraction of the true solution. As the cylinder accelerated, the reduced 
order variable scaled the mode to match the true solution. When the cylinder came to 
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(a) Analytical Solution (b) POD/ROM with 1 mode 

Figure A.4 POD/ROM/RG Contours of 4/ 

rest, the reduced order variable scaled to zero. The POD/ROM/RG with one mode was 
exact (Fig. A.4). The eigenmode could be captured at any time increment with a single 
snapshot. 

A.2.2 Deforming Grid Case. The three deforming grid sub-sub-cases were com¬ 
pared to determine if the type of deformation affected the accuracy of POD/ROM. 

Translating Sub-Case: In the translating sub-case, the cylinder was translated 
within the outer boundary to create compression and stretching of the cells. The primary 
mode of the translating sub-case was a scalar fraction of the analytical solution and ac¬ 
counted for 98% of the total eigenvalue energy. The other modes were radial in nature, 
due to the grid deformation maintaining straight radial grid lines (Fig. A.5). The remain¬ 
ing modes decreased in energy rapidly and were negligible by the 7th mode (10^^) (Fig. 
A.6). In the translating sub-case, POD/ROM only required four modes for a 99% accurate 
solution. Using fewer modes resulted in a slightly compressed solution (Fig. A.7). 

Rotating Sub-Case: In the rotating sub-case, the outer boundary was rotated 
with respect to the cylinder, creating skewness in the grid cells. As in the translating 
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(g) Mode 7 (h) Mode 8 (i) Mode 9 

Figure A.5 Translating Grid Sub-Case POD/ROM/DG Modes 
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Figure A.6 POD/ROM Modal Contribution (Rigid and Deforming Grid Cases) 





(a) Analytical Solution (b) POD/ROM With 2 modes (c) POD/ROM With 4 modes 

Figure A.7 Translating Grid Sub-Case POD/ROM/DG Contours of T 
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sub-case, the primary mode was a scalar fraction of the analytical solution and accounted 
for 99.7% of the total eigenvalue energy. The other modes were annular in nature, due to 
the annular sweep of the grid points as they rotated about the cylinder (Fig. A.8). The 
remaining modes decreased in energy even more quickly than the translating sub-case and 
was negligible by the 5th mode (10^®) (Fig. A.6). In the rotating sub-case, POD/ROM 
required only three modes for a 99% accurate solution. Using only two modes resulted in 
a small deviation in the areas of highest grid movement (the boundaries) (Fig. A.9). 

Translating/ROtating Sub-Case: In the translating/rotating sub-case, the grid 
was deformed by both translation within the outer boundary, as well as rotation of the 
outer boundary with respect to the cylinder. The translating/rotating sub-case was the 
most complex and required the most modes for an accurate solution. The first mode of 
the deforming grid POD/ROM had 98.7% of the total energy of the eigenvalues. The 
remaining modes decreased in energy quickly, but were still significant compared to those 
of the rigid grid POD/ROM (Fig. A.6). The first mode of the translating/rotating sub¬ 
case was similar to the rigidly translating case: a scalar fraction of the true solution. The 
remaining modes for the deforming grid POD/ROM were variations to account for the 
changes in the physical locations of the grid points (Fig. A. 10). Mapping the modes to an 
undeformed grid, the higher modes appear as corrections from the undeformed grid to the 
deformed grid. The translating/rotating sub-case POD/ROM required at least 7 modes 
for a 99% accurate solution. Using fewer modes resulted in a skewed solution (Fig. A. 11) 
where the principle mode is incorrectly mapped to the deformed grid. Using more modes 
resulted in a more accurate solution. 

If the POD/ROMs were run on a grid that differed from the training grid, the 
accuracy diminished. Two cases were tried, that of a synchronous clocking of translation 
and rotation, and that of an asynchronous clocking. If the boundary was started at zero 
degrees, the rotation was synchronous with the translation of the cylinder. Thus, the 
solution started at the undeformed/unrotated state and returned to that state at least 
twice more. However, when the solution started at the maximum rotation, the rotation 
was not synchronized with the cylinder movement and the solution never approached the 
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(a) Analytical Solution (b) POD/ROM With 2 modes (c) POD/ROM With 3 modes 

Figure A.9 Rotating Grid Sub-Case POD/ROM/DG Contours of \I/ 

undeformed/unrotated state within the time interval considered. In both cases, the grid 
never underwent the exact position of the training grids. 

In the synchronous sub-case, 11 snapshots (and hence 11 modes) could be used 
before the POD would become unstable due to linear dependance. The 10 grid modes of 
the system were very small compared to the initial fluid mode and had eigenvalues that 
were only 3% of the eigenvalue sum. By plotting the contribution of each mode to the 
maximum value of stream function, the relative magnitudes of the modes can be isolated 
(Fig. A. 12). 

In the asynchronous sub-case, 20 snapshots could be taken while maintaining linear 
independence. The additional number of snapshots (and therefore modes) was due to 
the more complex interaction between the asynchronous movement of the translation and 
rotation. The 19 grid modes of the system were very small in magnitude compared to the 
initial fluid mode, and had eigenvalues totaling only 5% of the eigenvalue sum (Fig. A. 13). 

In both sub-cases, the greater the number of modes included, the more accurate the 
solution was (Tbl. A.l). Even with the largest number of modes possible, the maximum 
error throughout the 100 time units was almost 10%. This error occurred approximately 
at midpoint in time between two snapshots, during the most rotated and deformed po- 
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Table A.l Maximum Percentage Error of Stream Function 


Synchronous sub-case 


Asychnronous sub-case 


# Modes 

% Max Error 

# Modes 

% Max Error 

6 

3380.50 

14 


8 

94.95 

16 

108.30 

10 

17.04 

18 

44.60 

11 

7.05 

20 

13.02 


Table A. 2 Effect on Error due to Mode Truncation 


Synchronous sub-case 


Asychnronous sub-case 


# Modes 

% Max Error 

# Modes 

% Max Error 

6 of 11 

1252.14 

14 of 20 

300.73 

8 of 11 

240.56 

16 of 20 

124.84 

10 of 11 

21.56 

18 of 20 

63.76 

11 of 11 

7.05 

20 of 20 

13.02 


sition. The average error over the entire grid of the stream function at this point was 
approximately 1.0%. 

Taking fewer than the maximum allowable number of modes (based on the limitation 
of the current computational method) resulted in a less accurate solution (Tbl. A.2). Mode 
truncation at the 99.9% level resulted in errors in excess of 1000%. 


A.2.3 Grid Density. Grid density was also examined. Solutions were calculated 
using the maximum number of modes for each sub-case using grid densities of / , J = 
10x10, 20x20, 30x30, and 40x40 (Tbl. A.3). As grid density increased, accuracy decreased. 
This was due to the deformation of each individual grid point being more severe. Each 
point had to move further relative to its original position in the more dense cases, in terms 
of index position. 


4 

rable A.3 Effect on Error due to Grid Density 


Synchronous sub-case 

Asychnronous sub-case 

IxJ 

% Max Error (11 modes) 

% Max Error (20 modes) 

10x10 

0.19 

0.32 

20x20 

0.60 

3.07 

30x30 

7.04 

13.02 

40x40 

22.47 

24.47 
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A. 3 Conclusions 


Overall, the use of POD/ROM with a deforming grid was reasonable. An algo¬ 
rithm must be created to ensure that accuracy is maintained, but the basic application 
of POD/ROM was sound. Of the above issues, the increase in the number of overall 
modes by the addition of distinct grid modes presents the main research area. As the 
number of modes increase, the solution time increases, decreasing the effectiveness of using 
POD/ROM. The additional grid modes cannot be truncated without significantly decreas¬ 
ing accuracy. The sample problem de-coupled fluid dynamics from grid dynamics. In the 
case of the rigid grid, the POD/ROM was able to return an exact solution with only one 
mode. In the case of the deforming grid, the POD/ROM was able to reproduce an accurate 
solution, but required more modes to do so. In all cases, the primary mode was similar 
regardless of whether or not there was deformation or not. This mode represented the fluid 
dynamics that was de-coupled from the grid dynamics. For the deforming grid cases, the 
remaining modes corrected the first mode to the deformed state of the grid. 

The magnitudes of the eigenvalues showed that POD/ROM was more sensitive to 
translation than rotation in this model problem. This was due to the greater scope of 
movement for all of the grid points in the translating case. The lower energy modes of 
translation and rotation were also significantly different. The grid modes tended to align 
themselves along the path of the grid motion. 

When translation and rotation were combined, the problem became even more com¬ 
plex and required more modes to solve accurately. The lower energy grid modes of the 
translating/rotating sub-case were not simple additions of the grid modes from the sepa¬ 
rate translating and rotation sub-cases. 

The grid density for any viscous problem will be an issue if a deforming grid is 
implemented. The grid density necessary for a viscous boundary layer is very large and 
any deformation will result in a large number of grid modes. As the number of modes 
necessary for an accurate solution increases with grid density, some method of minimizing 
the area of grid deformation will be necessary. One approach would be to limit the grid 
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deformation to a percentage change in relative position. If a grid point’s location does not 
change much in a single time increment, it could be kept static. 
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(a) Analytical Solution (b) POD/ROM With 4 modes (c) POD/ROM With 7 modes 

Figure A.11 Translating/Rotating Grid Sub-Case POD/ROM/DG Contours of 'll 



Time Untts 


Figure A. 12 Max value of Stream Function by mode 
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Figure A. 13 Max value of Stream Function by mode 
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Appendix B. Proper Orthogonal Decomposition Theory 

Also known as Karhunen-Loeve decomposition, principal components analysis, singular 
systems analysis, and singular value decomposition (1), Proper Orthogonal Decomposition 
(POD) generates a basis for the modal decomposition of functions and discrete data, 
such as experimental data. POD is one member of the class of representations known as 
orthogonal expansions (the Fourier series, or harmonic decomposition, is another example). 
The POD basis is linearly optimal in that it provides the most efficient way of capturing 
the dominant components of an infinite-dimensional process of finite dimension modes. 
The basis functions it yields are called empirical eigenfunctions, empirical basis functions 
and empirical orthogonal functions (1). 

The derivation of POD is recreated in a variety of works (see (1), (45), (14), (13)). 
However, the most complete derivation in the literature is presented in the dissertation of 
Newman (46), and it is from this source that the following is section is taken. 

B. 1 Derivation 

A stochastic process is a family (ensemble) of random variables or functions from 
some measurable space, from which a probability space can project into a state space. 
Suppose that is an ensemble of a stationary (no time varying mean values) stochastic 

process, defined on a spatial domain D. To obtain an accurate representation of the 
members of {X}, each is projected onto a set of candidate basis functions. Assume that 
X belongs to an inner product space (the linear, infinite-dimensional Hilbert space, £2 (D), 
of square integrable (complex-valued) functions), the orthogonal representation of 

00 

X {t,x) = ^anit)(l)^{x), (B.l) 

n 

where is the orthonormal basis and {on} the uncorrelated set of random variables, 
in time. A reduced-order model (ROM) can be obtained by truncating the series equation 
(B.l) to N terms. 
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The assumptions of the inner product space affect how POD is developed. A Hilbert 
space 7-f is a vector space over M or C together with an inner product (•, •) which is complete 
as a metric space. The norm is defined as ||0|| = {(j), (j)) for 0 G W and the metric is 

defined as = ||0 —-011 for 0,-0 G 7-f. A set of vectors 0 in a Hilbert space Ti is 

an orthonormal set if any two distinct vectors 0 i ,02 ^ ® ^^6 orthogonal, i.e., (0,0) = 0, 
and in addition, ||0|| = 1 for each 0 G 0. The £2 [a, &] space is that of complex-valued 
(or real-valued) Lebesgue-measurable, square integrable functions on a domain D, and is 
generally an infinite-dimensional Hilbert space with an inner product on £2 {D) of 

U^9)c2{d) = j f{x)9{x)dx. (B.2) 

Newman provides a detailed derivation of the Hilbert space, TL, and its properties 
as they relate to POD. The development of POD requires an orthonormal set of basis 
functions, {0}, for the £2 {D) subspace. Because it is a Hilbert space, there are an infinite 
number of orthonormal basis sets. However, Sirovich (47) provides a method of selecting 
a particular choice. 

Theorenn 1 Sirovich: //{0} is a complete orthonormal set of functions in £ 2 ( 7 ?) such 
that the expansion can be written as equation (B.l) where both of the sets {an} and {4>n} 
are orthonormal, i.e., 

iflki ^l)C2{D) ^ (^) (^)] l®'^) 

and 

{(t>k, ^l)c2iD) ^ 

with E[-] as a probabilistic averaging operation (1), then {(fn} found in a spectral 
decomposition of the two-point spatial correlation function 

R{x,y) = '^<Pn{x)(l)n{y) ■ (B.5) 

n 

This then provides an attractive candidate for the orthonormal basis to use in the 
expansion, and the goal becomes to determine if such a basis exists. The existence of the 
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expansion, equation (B.l), that satisfies the above conditions is guaranteed under certain 
conditions by the Karhunen-Loeve expansion theorem (presented below). Moreover, the 
Karhunen-Loeve expansion theorem proves a method for constructing the orthonormal set 
of basis functions. The proof depends on the existence of the spectral decomposition of 
the spatial correlation function, which is guaranteed under certain conditions by Mercer’s 
theorem. 


Theorenn 2 Mercer: Let be a continuous, Hermitian symmetric, nonnegative def¬ 

inite function on [a, 6] x [a, 6]. If {(j)^} and {A^} are a basic system of eigenvectors and 
eigenvalues of the integral operator with kernel k{-,-) then ys,t £ [a, b], 

OO 

k{s,t) = ■ (B.6) 

n=l 

The series converges absolutely and uniformly on [a, 6] x [o, b]. 

With the Mercer theorem, the Karhunen-Loeve theorem defines the continuum of 
random variables by a countable number of orthonormal random variables. 


Theorenn 3 Karhunen-Loeve Expansion: Let {Xt,t G [a,b]} be a zero-mean, quadratic- 
mean continuous second-order stochastic process with covariance function R{t,s). The 
process Xt has an orthogonal decomposition 

N 

X{(jj,t)= lim ^/xlai {(Jj) (j)i (t) te[a,b], (B.7) 

i^oo 

with an associated averaging operation 


E ttj] — 


(B.8) 


and 


fb 

‘^j)c2[a,b] ~ J ~ 


(B.9) 
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if and only if the {(j)i,(j) 2 , ■■■} o,re the orthonormal eigenfunctions and the {Ai,A 2 ,...} are 
the eorresponding eigenvalues of the integral operator with the kernel R{-, •) i.e., 

R{t,s) (j)^{s) ds = Xi(j)^{t) t G [a, 6] i = l,2,... (B.IO) 

In that case, the series converges uniformly on [a, b]. 

The coefficient functions {</>„} that correspond to the non-zero eigenvalues {An} are 
called the empirical eigenfunctions. They form an orthonormal basis for the subspace 
C 2 {D). The uncorrelated random variables {an} are given by 

On (w) = (^\/{An}^ J (j)n{t)X{‘^,t)dt, (B.H) 

and form the desired orthonormal basis for R. 

To summarize the theorems, if there is a stochastic ensemble and if it can be expressed 
as equation (B.l), then a ROM can be created by truncating the orthonormal basis vectors 
{<fn}- B the space is a Hilbert space and the inner product is defined as equation (B.2), 
then there are infinite orthogonal basis sets. Sirovich’s Theorem says that one of those 
basis sets can be defined using a probabilistic averaging function and the spatial correlation 
function (kernel). Mercer’s theorem defines a kernel and Karhunen-Loeve expansion finally 
defines the orthonormal basis in terms of the original stochastic ensemble and the kernel. 

B.2 Discretized Sample Spaces 

In discrete applications, the ensemble is sampled in time, space, or both. The dis¬ 
cretized Karhunen-Loeve expansion relies on the spectral theorem for real symmetric matri¬ 
ces which is the finite-dimension analogue to Mercer’s theorem. If the set, Xi (u) ,i ^ [a,b], 
is a zero-mean scaler-valued, discrete-parameter, second-order process with covariance ma¬ 
trix R, then the matrix R can be expressed as R {j, k) = E [XjXj]. The matrix R is real, 
symmetric, and non-negative definite. The spectral theorem states that every matrix R 
can be diagonalized by an orthogonal matrix. That is, there exists an orthogonal matrix 
and a real diagonal matrix A such that R = The spectral decomposition of R 
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can be expressed as 

OO 

R = Y.^n^n4'l, (B.12) 

n=l 

where each vector is the n-th column of 

T heorenn 4 Sampled Data Karhunen-Loeve Expansion: Let {Xj (a;), i G [a, 6]} be a zero- 
mean scalar-valued diserete-parameter, second-order process with covariance matrix R. The 
process Xj has an orthogonal deeomposition given by 

N 

X {w, i) = lim ^ \/Kan (a;) k = 0,l,... (B.13) 

N—^oo 

n=\ 

with 

E[aiaj]=6ij, (B.14) 

and 

(t>j) = = Sij, (B.15) 

if and only if {fn} the orthonormal eigenveetors and {An} are the corresponding eigen¬ 
values of the matrix R: 

Rfn = K4>n- (B.16) 

B. 3 Optimality 

The optimization of POD can be expressed as a minimization of error or, equivalently, 
as the maximization of the projection of any member of the stochastic ensemble onto 
the subspaces spanned by the most energetic members of the empirical basis. Consider 
the ensemble {X (cv,-) : [o, 6] —*• G 0} and let {4>n} be the empirical eigenfunctions 
with corresponding eigenvalues, {An}. Let {X (w, t) ,uj & D,t & [a, 6]} be a member of the 
ensemble with POD 

N 

X(c^,t) = J]6n(ch)0nW, (B.17) 
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where bn (lu) = \/X^an (a;). Let {V^^} be an arbitrary orthonormal set such that for some 
random variables {cn} 


N 

^ t) = ^Cn (w) Ipn (t) ■ 


(B.18) 


n=l 


Then, for each truncation index N 


N 


N 


n=l n=l 


N 

n=l 


(B.19) 


and equivalently 


N 

n=l 


\\bn 


N N 

E^">E® 

n=l n=l 



(B.20) 


To examine the minimization of error, on average, between the members of the 
ensemble and the truncated orthogonal expansion, observe that the minimum of the mean- 
square error 



00 

2 " 



N 



Xt-'^Cn (OJ) i)n {t) 

n=l 


= E 


+E^ 

n=l 

\\Cn{u))i)^ (t)f 


is achieved when '^n=i E [(cn (<^) i^n (^))] is maximized. 


N 

-2 [(Cn (w) f/’n (t))], 

n=l 


(B.21) 


B .4 Method of Snapshots 

From a practical point of view, data is sampled at discrete values of time. The 
calculation of the spatial covariance matrix R = E [XjX/"] must use discrete samples or 
snapshots of data. These sets of snapshots form a snapshot matrix, S, where the sampled 
ensemble data, {X (ti ), X (^ 2 ) , •••}, is placed into columns. If the process is assumed to 
be ergodic or stationary (time average equals ensemble average for each fixed value of X) 
then the averaging function can be expressed as 


R = E [XtXf 


lim 

T—^00 


1 

M 


XtXf dt. 


(B.22) 


B-6 



R can be further approximated as 


R^R= lim ^ E X {tn) X {tnf . (B.23) 

M—^oo IVi 

n=l 

Since there are only a finite number of M snapshots available, the approximation, R, is 
further reduced to 

1 “ 

= ^ (B.24) 

n=l 

Rm can be expressed in matrix form as 


Rm = (B.25) 

where S' is an x M matrix containing M snapshots of X (t), with the data stored 
as column vectors in the matrix. It is common for authors to ignore the ^ factor and 
absorb it into the generation of the eigenvalues in the solution of the eigen-problem of 
equation (B.16), i?M0n = ^n4>n- This particular problem can be solved by Singular Value 
Decomposition, yet it requires the calculation of the N x N matrix Rm to be explicitly 
carried out. 

In fact, the Method of Snapshots does not use the spatial covariance matrix, Rm, 
but instead a temporal covariance matrix, Cm- Assume that the snapshots are linearly 
independent vectors, i.e., the data matrix S has full column rank, and define the temporal 
covariance M x M matrix 

Cm = j^S^S. (B.26) 

Starting with a similar eigen-problem as equation (B.16), and 

CM'4^n 

and using the definition of Cm and eigenvectors and values of {V’n} and 

i G M, (B.27) 
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S fJ>n'4^n 


(B.28) 


The eigen-problems are equivalent if and only if 

K = Stpn i e M. (B.29) 

Resulting in, 

= Rm (pn (B.30) 

The eigenvalues {/r„} for Cm correspond exactly to the non-zero eigenvalues {An} of 
Rmi but the computational requirement for Cm is much smaller, as for any practical prob¬ 
lem M « N. The relationship of the empirical eigenvectors {Cn\ the eigenvectors 
{Cn} of Cm are given as 
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Appendix C. Full-Potential Equation 

The derivation in this appendix is an expansion of the work presented by Shankar, et al. 
(33). This expansion shows a full step-by-step development of the discretized equations 
used in this method. The full-potential equation is an approximation to the Euler equations 
for which the flow is assumed to be irrotational, free of entropy production (36). The 
potential equation has one fluid variable, T, defined as. 


d'S 

dx 

dy 


(C.l) 

(C.2) 


where u and v are the fluid velocity components (Note that (j) is typically used as the 
potential variable but has already been assigned in this document). In a 2-D, body-fitted 
coordinate system represented by r = t, ^ ^ and 77 = ri{x,y,t), the unsteady 

full-potential equation is written as 



r 


+ 




(C.3) 


where the (R)^ refers to the derivative of R with respect to U and V are the contravari- 
ent velocities, J is the Jacobian of the metrics (j = ~ ^yVx) the density, p, is 

represented as 


P = 


= 1 


7-1 


Ml [2T, + {U + U + (^ + Vr) - 1] 


1 

7-1 


(C.4) 


The contravarient velocities are defined as follows: 


U = -|-J-0124/,,, 

V = T/r + 

The transformation grid metrics an, an, and an are defined as: 

On = ^x + ^y 


(C.5) 

(C. 6 ) 


(C.7) 
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®22 11X Vyi 


®12 ^xVx 


(C.8) 

(C.9) 


The equations above can be solved by a variety of schemes. The scheme selected was 
developed by Shankar, et, ah (33) and has the benefits of being very stable and easy to 
implement with POD. It uses an implicit Newton’s method to solve Eqn. C.3, in discrete 
form, for a robust and efficient numerical solver with effective treatment of boundary 
conditions. In summary, the scheme is as follows: Eqn. C.3 is represented as 




(C.IO) 


where T is the unknown to be solved for at every grid point in the domain for time step 
(n + 1). The Newton iteration for the solution to Eqn. C.IO is, 


E(T,) + 


dF 


= 0 , 


(C.ll) 




where is the current guess for T at the {n + I) time step. At convergence (typically 
less than 10^®), AT = T — will approach zero. 

The time derivative, and spatial derivatives, and ; must be treated 

in a manner consistent with the Newton iteration of Eqn. C.ll. The method assumes 
that any variable may be expressed as a deviation from an assumed state, f = f (T^.) + 
AT 


C. 1 Treatment of 

A first-order time discretization is applied, 


m" - (5) 






At 


(C.12) 


where At is the change in time from step (n) to (n-l-l), and (i,j) are the grid-point 
coordinates. The unknown in Eqn. C.12 is and can be evaluated in a time-linearized 
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manner similar to Eqn. C.ll: 




(C.13) 


The differential operator, ( ^ ) , can be derived as follows: 


dp 


d-^ 


dp 


^ 1 AM/. = ^AM/. + ^AT^ + 






dp 


dp 


(C.14) 


Assuming the grid metrics, and, ry. are small, then: 



[24-. + (U + 5.) S-f + (V + 77.) S', 
[247. + (a + 5.) 'Pf + (V + 7,^) S', 
[2S'. + (U + 5 .) »£ + (V + 7,^) S', 



(C.15) 

(C.16) 

(C.17) 


This results in: 


dp 

dp 

dp 

The speed of sound can be defined as, a 



(C.18) 

P Ml, 

(C.19) 

r OO TT 

2 ’ 

p Ml, 

(C.20) 

p7-i 2 


resulting in: 


dp 

dp 

dp 





(C.21) 

(C.22) 

(C.23) 
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The time derivative is approximated by and eliminating the constant, 2, results in a 
final form that is evaluated at T = for the p, J, U, and V terms: 


+ 




At 


drj 


AT. 


(C.24) 


J '1'=^* 


It may be differentiated using second-order, central-space operators on the half grid point: 




/ 






‘■('f.)y 


A^i 


At 






(C.25) 


Finally the time derivative is determined to be: 


\ J ) T Xt V J / i,j \ J) i,j \ JA-r^a j 


/ 

V 




-^+ 
['A'P. i, ,-A'P, 


-l.j) 


} 


(C.26) 


This is reduced to: 



( P{^*) \ 

\JATa{^^f ) 


f JAT^a{9^f\ 



V Jij 

[V 

[jJiJ 


+ 


AT 


- U (T=, 




2 


V{^.\,At 


. (A^i.,+i-A^i,,,i) 
2 


(C.27) 


The factor, () that appears in Eqn. C.27 is used in the later development of 
the Approximate Factorization. 


C.2 Treatment of 


Again, using the Newton procedure from Eqn. C.ll, the spatial derivative is 

factored as: 


pU 


J 


d 


— = — f H—^AT 


df 


where: 






9T 




(C.28) 


(C.29) 
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By chain rule, 


(C.30) 


^ = -( r r 9p \ 

J \d^ ) ’ 


where is given by Eqn. C.15 and is developed as follows: 


dU 

dU 


dU 








ni 


V 


aiiAil/^ + ai2A\l'^. 


(C.31) 

(C.32) 


For mixed flows, where elliptical and hyperbolic regions coexist, the inclusion of leads 
to a pentadiagonal matrix. To preserve a tridiagonal form for efficient matrix inversion, 
the term U ^ appearing in Eqn. C.30 is neglected. This is acceptable in that U^ is 
multiplied by AT which approaches zero. This results in the spatial derivative term, 


J )i 


d 

di 




J V^T ) 


(C.33) 


This expands to: 


_ d fp(T, 
j)^ diX J 


, 9AT 5AT 
U (T=^) + -h 012 


di 


drj 


(C.34) 


Substituting the definition of U (T^,) yields: 

+ on (T=, + AT)^ + 012 + AT)^ 


(pU\ 

d 

f p(^*) r 

\j), 


1 ^ ^ 


(C.35) 


The spatial terms are differences using second-order, central difference operators on the 
half grid-points, 

/ nTT\ /nTT\^+^ /nTT\^+^ 

(C.36) 


n 




p£ 

J 




Numerically, the flux terms are calculated in a single sweep and the i + ^,j points. After 
substitution, the spatial derivative becomes: 
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pV \ "+1 



[(’I'.ij + l - ’I'.i.j) + (A'l'ij - 

H- 4- - — + + l + 1 — 

+-(A'I'i+1,3 + 1 + - AS-i.ij + i - Avl/i^i_j) 



(C.39) 


C.4 Summary of Solution (j)^ + + (^) ~ 


The overall solution is obtained by substituting the above equations in the Newton 
iterative equation, 

(A^) = 0. (C.40) 



(C.41) 
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F(9,) = 


where /3 = 




-1,3 + 1 
1-1,3 
-1,3-1 
1 


JATa(^«)^ 


( pCH, ) ) 

{ £12 ^ 

('Pd'.)'! , 

V .1 ) 

V 4 

1 J J ' 

^ 14 + 1,3 

+ {^){ 

^ ^..3+1 

/'p(+.)'l 



V J J 

i ^ b+i. 

j \ J ) 

al2 1 

(Pi'S.) 


1 h+1 

,3 V 1 

K 4 


I 'I 

_ I'A+ii'i r (°i2) i 




(“ii). 


+ 3,3 V J 
+ 




(“ll),^l ,- + (^1 (“22) 


1 1 / 1+1,3 




1 


. J j 


+ 

3 

' (^ 12 ) 


j 




J—y (^ 22 ). 




(‘^ 11 )+ 




(^)(^)._i., + (^)(^) 

pOtii) _ Z'iL'l" 1 




'pi+ii) 1 ( 212 !'' 
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C. 5 Approximate Factorization 

The Jacobian matrix, (§^), is large {N x N) but sparse and can be simplified using 


an approximate factorization creating two linear operators, 


Lr = 


Lri — 


1 + ArU^ + 

1 + AtV-^ + 
op 


JAr^a(T^.)^\ d /' p d 



d f p d 
I dp dp 


(C.43) 

(C.44) 


where and are {N x N) tridiagonal matrices that are solved with successive ^ and 
p sweeps 


L^L^AT = R, 

(C.45) 

L^AT = R, 

(C.46) 

L^AT = AT, 

(C.47) 


where AT is an intermediate solution and R = F (T^.) from Eqn. C.ll. 

The linear operators are discretized in a fashion similar to the full development with 
second-order, central difference operators on the half grid points: 
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= 


A'i/ + AtU — + 

9? 


JAr^a 

p('i'.) 


d p d p d‘^ 

— ail-A^' + ail-A^' 

di j ae j di^ 




AS'i.i + ^ (A'I'i+1,2 - Avl/i_i,p 


L^A'^ = 

( JAT2a(<I-.)2\ 1 




+ V ) 1 

1 +(“ll5),_7A'I'i+l,2 J 



(C.48) 


(C.49) 


A4' + AtV — A<^ + 
dt] 


JAT^a{<if^Y 

p(9,) 


a p a p a^ 

— 022 -—A* + 022-;^A'I' 
or] J orj J dr]‘^ 




i,3 + l - A'l'i^j^i) 


( JAT2a(<I-,)2\ 1 




^ V ) \ 

1 +(“""j).,2 

(A'l'i.j + i - 2AvI>i_j + J 



(C.50) 


(C.51) 


C.6 Boundary Conditions 

The boundary conditions are broken into three groups: far-held, body surface, and 
unsteady wake . 


C.6.1 Far-Field. The far-held boundary condition uses Riemann invariants that 
correspond to the positive characteristics with respect to the inward normal. Thus, for the 
inhow, outhow and far held boundary the positive characteristics are: 


U 2 

, — H- 

YaE 7-1 

U 2 

7-1 

V 2 

7-1 

These can be differentiated with respect to T: 


const, 

const, 

const. 


I dU 2 da 

VoITM + 7-1^ 
I dU 2 da 

7^^ + 7-1^ 
1 dV 2 da 

^/a^ 9T 7 — 1 5T 


(C.52) 

(C.53) 

(C.54) 


(C.55) 

(C.56) 

(C.57) 
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Substituting the derivatives for and and setting the residual to the difference in 
the freestream and the calculated value at the boundary yields, 


7^ (aii^ + 012^) AM/ + ^ + U-§^ + am;' 

“ ( 7 ^ freestream ~ 

“12 am/ + ^ + 7 ^ + V-^'^ AM/ 

“ (“7^ + ~“) freestream ~ (“7^ 

(“12^ + “22^^) A^f + ^ + 7^ + V-§pi^ 

~ (“7^ 7d:“) freestream ~ i~7^ 7^:“) 


(C.58) 


(C.59) 


(C.60) 


The scheme is discretized using first-order forward or backward differences for the normal 
derivatives and second-order central differences for the parallel derivatives, 

Inflow / Outflow: 


(-y —l)A-r 



+ V 




(C.61) 


Far Field: 


(f -1) At 




-u) 

^ + ■. 


nl 


: + V 


^ freestr 


(C.62) 


C.6.2 Body. The body boundary condition requires flow tangency for inviscid 
flow. The contravarient velocity, V, must be zero at the body. By definition. 


By assuming: 


F 


rjr + 012^'^ + 022 ^'r? = 0, 


Vt 

022 


012 

022 





(C.63) 

(C.64) 


(C.65) 
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where the i,j — ^ point is a ghost point within the body, the governing equation can be 
modified for a body point. For the Approximate Factorization scheme, the implementation 
of the body boundary condition must be done in both of the linear operators, as 

well as in the right hand side term, R. The linear operators are modified as follows; 


Li 


L 




1 + A,t/| + 

\ pi^,) 


0-22 


an - 




aji 

022 


p d 


(C.66) 

(C.67) 


They are then discretized as above. The right hand side term, R, is modified in the 
term. 


pU 

J 





d^ 1 

1 J 



d^ 1 

[ J 



d^ 1 

1 J 


U ('1'=^) + Oil—-h 012 


U (IP*) + Oil 


9AT 


dp 


di 


Vt 012 


+ 012- 


dA-^\ 


022 022 


d^ J 


oii-^ 


ah \ 


022 


022 / d^ 


The 



term cancels. 


(C.68) 

(C.69) 

(C.70) 


C.6.S Unsteady Wake. The wake cut behind the trailing edge of the airfoil must 
be properly modeled in an unsteady flow. The modeling is done through the vorticity 
convection equation. 


Ft + [oii'S^^u - ^t,i) Li,u + {^t,u + Ct,/) '^i,u — 0 (C.71) 

where the u, I subscripts indicate upper and lower wake boundaries. The vorticity con¬ 
vection equation assumes that ~ 0, where [/] is the jump in the quantity / across 

the wake. Also, 012 is assumed to be zero for all wake points by proper construction of 
the wake grid. The vorticity equation is integrated via a direct marching scheme from the 
trailing edge to the outflow boundary condition. It is begun with the assumption that 
F = ~ at the trailing edge. Finally, after the vorticity distribution is determined. 
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the lower wake potential function values are overwritten, 




(C.72) 


To maintain symmetry and preserve accuracy, on every odd time count the potential 
function on the lower wake cut is preserved and the upper overwritten. 


Tt + = 0 . 


(C.73) 


In addition to the vorticity distribution, the unsteady solution also requires on the 
wake points. This is determined via a Taylor’s expansion, 




+ T + 




(C.74) 


where: 


[^vv] ~ 


paiiTj \ 

V J J 


(C.75) 


C.7 Density Biasing 

All of the density terms in and are density biased to account for shocks. 

Several density biasing techniques were proposed by Shankar, et. al. The directional flux 
biasing was selected to provide a fully rotated solution. 

P= ^ , (C.76) 

where q is the local total velocity and is defined to be (pq)^ = PQ — P*Q* H q > q* and 
(pq)^ = 0 if q < q*. For U > 0, the positive sign and forward differencing is used. For 
U < 0 the negative and backwards differencing is used. The quantities p*q* represent sonic 
values of density and total velocity. These are given using the speed of sound relationships 
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as: 


(q*? = 


1 + [1 - 2V1/, - - 2r^,%] 


1+1M'^ 

2 -^*^^00 


P = 




(C.77) 

(C.78) 
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Appendix D. Validation 

Validation of the Shankar full-potential code used for the full-order solver was done by 
comparison to experimental surface pressure data and the results of previously validated 
computer codes for two cases; steady flow over a 2-D stationary bump (40) and unsteady 
flow past a NACA 0012 airfoil (48). For additional comparison, inviscid and viscous results 
from the COBALT computer code (49) were compared to the steady case. 

D.l Steady Flow 

The results of the Shankar code compared well with the experimental subsonic data 
and with the subsonic and transonic COBALT code. In the transonic cases, the Shankar 
code provided accurate integrated surface pressure as compared to the experimental data. 
The experimental data was taken for a smooth bump in the Langley rectangular high¬ 
speed wind tunnel (1952). Pressure measurements were taken for two bumps of thickness 
to chord ratios of 0.05 and 0.15 (Note that Lindsey and Daley list ratios of 0.1 and 0.3, 
but their definition of chord length is from the midpoint to leading edge). Mach numbers 
from 0.25 to choking flow were tested at zero angle of attack. The tunnel had a four by 
eighteen inch test section, with four-inch models. The tests were made on one surface 
of each model. The models were placed in the section together, back-to-back, assuring a 

zero angle of attack and dividing the channel into two separate test sections. The upper 

boundary was approximately 2.25 chord-lengths from the surfaces. Each model had one 
chord-wise row of static pressure orifices installed on the model surface at the semi-span 
station. The profile of the bump was made smooth such that reasonable agreement could 
be made between experimental data and potential flow at low speeds. The equations of 
the X and y coordinates of the surface were (50), 

X = cos 6* — ^ (cos0 — cos 30), (D.l) 

y = ^ (3 sin 0 — sin 30), (D.2) 

where 0 runs from 0 — vr and t is the thickness at the midpoint of the bump. 
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The results from the 0.05 ratio model were used for comparison to the Shankar 
code. The results of the Shankar code were slightly faster in maximum Mach numbers 
than the experimental data for the subsonic flow. This was to be expected, as the full- 
potential equation does not take into account the effects of viscosity. The increased velocity 
translated into a decrease in pressure on the surface of the bump (i.e. larger negative Cp). 
The same result was seen in the COBALT (Euler) cases and the two computational codes 
were within 5% of each other (see Fig. D.l). The fully subsonic cases show good agreement 
between potential, Euler, and experimental results. The subsonic profiles were symmetric 
for the Euler and potential codes. The experiential profile was not symmetric due to viscous 
effects over the rear half of the bump. At an inlet Mach of 0.759, the flow becomes transonic 
at the top of the bump. Both the Euler and the potential results begin to show a shock, 
however the experimental result does not. In the experimental results, the shockwave does 
not show up in pressure data until about Mach 0.789. This is due to the effects of the 
boundary layer altering the effective shape of the bump. At Mach numbers greater than 
about 0.8, the experimental flow fully separates behind peak of the bump. Therefore, it is 
impossible to compare the potential results to the experimental data. However, the Euler 
and potential results continue to be very similar (Fig. D.2). 

D.2 Unsteady Flow 

The results of the Shankar code compared well unsteady experimental NACA 0012 
airfoil data. The experimental data was transonic with a weak shockwave on the body of 
the airfoil. Oscillation was forced in both pitch and plunge and pressure data was recorded, 
then converted to non-dimensional coefficient of pressure and lift coefficient. The Shankar 
code showed excellent comparison to the experimental data, in fully developed flow (Fig. 
D.3), and matched the original data presented by Shankar et. al. (33). The Shankar 
solver clearly reproduced the Euler results and reproduced the experiential results within 
the limitations of the full-potential equation. 
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(a) Mach = 0.226 (b) Mach = 0.519 



PfmrlCtv^^ FmiCtv^^ 


(c) Mach = 0.680 (d) Mach = 0.759 

Figure D.l Sub-Sonic Cases 



PmrlCtTTltU ' PmrlCtTTl^ 

(a) Mach = 0.789 (b) Mach = 0.829 

Figure D.2 Trans-Sonic Cases 
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(a) Lift Coefficent (b) Moment Coefficient 

Figure D.3 Unsteady Case 

Table D.l PAPA Structural Model Parameters For Validation 


Xcg 

Xa 

II 

TUry 

r2 

' a 

t^s 

Moo 

ao 

0.5 

-0.25 

0.0 

0.2 

0.25 

125 

0.8 

0 


D.3 Structural Model 

To validate the loosely coupled structural and fluid model, flutter boundaries were 
evaluated for the NACA 0012 airfoil. Data from the ENS3DAE CFD code and the TVD- 
ntiAE CFD code were used as baselines for comparison (51). The flutter boundary was 
expected between u = 6.6 and u = 7.0. The structural parameters selected were based on 
those used in the references (Tbl. D.l). 

Time integration was used to bracket the flutter onset speed. The results showed 
reasonable comparison to the TVDntiAE and ENS3DAE codes. The flutters onset speed 
was found to be between u = 6.8 and u = 6.9 (Tbl. D.2). The flutter onset stabilized after 
several hundred non-dimensional time units (see Fig. D.4). 
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Table D.2 


Flutter Onset Cases 


u 

Comments 

6.0 

Damped 

6.5 

Damped 

7.0 

LCO 

6.75 

Damped 

6.9 

LCO 

6.8 

Damped 



Figure D.4 Flutter Onset 
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Appendix E. Pitching and Plunging Airfoil Mode Shapes 

The following are mode shapes for the individual POD/ROM of the Pitching and Plunging 
Airfoil application. They are presented for greater detail in examining the nature of the 
mode shapes when varying pitch, plunge or maximum angle of attack. The first four modes 
are presented for several different input parameters: Figs. E.l to E.5 A = 0, a = 1.0 deg, 
/ = 0.02 — 0.1, and 20 modes; Figs. E.6 to E.8 A = 0 — 20, a = 1.0 deg, / = 0.06, and 20 
modes; Figs. E.9 to E.ll A = 10, a = 0.0 — 2.0 deg, / = 0.06, and 20 modes. 
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Figure E.l 


Individual POD/ROM Modes {h = 0, a = 1.0 deg, / = 0.02, 20 modes) 






Figure E.2 Individual POD/ROM Modes (h = 0, a = 1.0 deg, / = 0.04, 20 modes) 





















Figure E.7 Individual POD/ROM Modes {h = 10, a = 1.0 deg, / = 0.06, 20 modes) 



Figure E.8 Individual POD/ROM Modes {h = 20, a = 1.0 deg, / = 0.06, 20 modes) 



Figure E.9 Individual POD/ROM Modes (/i = 10, a = 0 deg, / = 0.06, 20 modes) 



Figure E.IO Individual POD/ROM Modes {h = 10, a = 1.0 deg, / = 0.06, 20 modes) 




Appendix F. Summary of Computational Runs 

The location and nomenclature of the archived computational runs are described in this ap¬ 
pendix. Archives are presented for each of the three model problems: Translating Cylinder, 
Oscillating Panel, and Pitching and Plunging Airfoil. 

The Oscillating Panel zipped archives are named based on the amplitude of deflection 
and the percentage of deforming area, i.e. DG075100 would be a maximum amplitude 
A = 0.75 and a 100% deforming domain. Within each archive the files are broken down in 
the following fashion: filename_f_.plt. The frequency of oscillation, /, is used to differentiate 
the separate files. Restart files are not used for the Oscillating Panel problem. The 
following zipped archives are stored for the Oscillating Panel problem: 

• DG05100.zip 

• DG075100.zip 

• DG10100.zip 

• DG1075.zip 

• DG125100.zip 

• DG15100.zip 

• DG2100.zip 

• DG1050.zip 

To execute a particular run, the user must locate the appropriate amplitude and 
deforming area. Within the zipped archive, the snapshot_/_.plt and the gridsnap_/_.plt 
files must be copied into the directory with the executable. Blended runs must have their 
snapshot files recreated using option (3) of the executable. Multi-POD runs must have all 
blended snapshot files stored in the same directory and then identified during execution. 

The PAPA zipped archives are named based on the amplitude of deflection, i.e. AlO 
would be a maximum amplitude A = 0.1. Within each archive the files are broken down 
in the following fashion: filename./_a_.plt. The frequency of oscillation, /, and AoA, a, is 
used to differentiate the separate files. Restart files are named in the following convention; 
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restart_M_o_.plt, where M identifies the inlet Mach number and a the initial AoA of the 
airfoil. All airfoils are assumed to start at a plunge depth of /i = 0. The following zipped 
archives are stored for the PAPA problem: 

• A005.zip 

• A010.zip 

• A015.zip 

• A020.zip 

• RestartM0510.plt 

• RestartMOSO.plt 

• RestartMOSOl.plt 

• RestartMOSlO.plt 

To execute a particular run, the user must locate the appropriate amplitude and 
deforming area. Within the zipped archive, the snapshot./_a_.plt and the gridsnap_/_o_.plt 
files must be copied into the directory with the executable. In addition the appropriate 
restart file must be copied to the same directory. Blended runs must have their snapshot 
files recreated using option (3) of the executable. Multi-POD runs must have all blended 
snapshot files stored in the same directory and then identified during execution. 

The source codes are stored in separate directories, with each code variation the 
name of a different sub-directory. The files types are stored within the zipped archive are 
identified in Tbl. F.l. 

An example case is stored in a separate directory, including all necessary data files 
and executable code. The example is a forced pitching and plunging airfoil at an oscillation 
frequency of 0.1, a plunge depth of 0.1, and a pitch angle of 2.0 deg. To execute the example, 
perform the following steps: 

1. Run Unsteady.exe 

2. Select option 5 
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Table F.l Archived File Names and Descriptions 


Eile Name 

Description 

fulljc_.plt 

full: Full-order results 

cljx;_.plt 

cl: Full-order lift (normal) coefficient values 

reduced_x;_.plt 

reduced: POD/ROM Results 

clredjc-.plt 

clred: POD/ROM lift (normal) coefiienct values 

phifield_x;_.plt 

phifield: Modal matrix, plotted on the computational domain 

snapshot _x_plt 

snapshot: Snapshots for individual amplitude and frequency combinations 

gridsnap_x_.plt 

gridsnap: Grid snapshots for individual amplitude and freuqency combinations 

error_x;_.plt 

error: Time accurate error between full-order and POD/ROM resutls 

maxer r or jx;_. pit 

maxerror: Peak error values over time 


(a) Select sub-option 1 

(b) Enter initial plunge depth of 0.0 

(c) Enter initial AoA of 0.016 

(d) Enter restart file name of RestartM05A0016.plt 

3. Select option 6 

4. Select option 4 

(a) Select sub-option 1 

(b) Select sub-option 1 

(c) Enter snapshot file name of snapshot.pit 

(d) Enter number of modes 25 

(e) Select sub-option 1 

(f) Enter number of modes 25 

(g) Select sup-option 1 

5. Select option 7 

(a) Enter number of iterations 10000 

(b) Enter number of oscillations 10 

(c) Enter plunge depth 0.1 
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(d) Enter pitch angle 2.0 

(e) Enter delta time 0.01 

(f) Enter number of sub-iterations 5 

(g) Enter accuracy of sub-iterations 0.000001 

(h) Enter data recording spacing 500 
6. Select option 99 for program end 
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